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Soon-Jo  Chung,  Sc.D.,  University  of  Illinois  at  Ur bana- Champaign 
Program  Monitor:  Dr.  Willard  Larkin 

Executive  Summary 

There  is  a  considerable  interest  in  developing  robotic  aircraft,  inspired  by  birds,  for  a  variety  of  mis¬ 
sions  covering  reconnaissance  and  surveillance.  Flapping  wing  micro  aerial  vehicle  (MAV)  concepts 
have  been  put  forth  in  light  of  the  efficiency  of  flapping  flight  at  small  scales.  These  aircraft  are 
naturally  equipped  with  the  ability  to  rotate  their  wings  about  the  root,  a  form  of  wing  articulation. 
This  AFOSR  Young  Investigator  Program  (YIP)  project  covers  comprehensive  research  problems 
concerning  the  performance,  stability  and  control  of  tailless  robotic  aircraft  with  articulated  wings 
in  flapping  and  gliding  flight. 

Neurobiologically  Inspired  Control  of  Flapping  Flight.  Significant  progress  has  been 
made  in  developing  a  rigorous  mathematical  model  of  a  neurobiologically  inspired  in  the  form  of 
central  pattern  generators  to  control  flapping- flight  dynamics.  A  rigorous  mathematical  and  control 
theoretic  framework  to  design  complex  three-dimensional  wing  motions  has  been  developed  based 
on  phase  synchronization  of  nonlinear  oscillators.  In  particular,  this  work  shows  that  flapping-flying 
dynamics  without  a  tail  or  traditional  aerodynamic  control  surfaces  can  be  effectively  controlled  by  a 
reduced  set  of  central  pattern  generator  parameters  that  generate  phase-synchronized  or  symmetry¬ 
breaking  oscillatory  motions  of  two  main  wings.  Furthermore,  by  using  Hopf  bifurcation,  this  work 
shows  that  tailless  aircraft  alternating  between  flapping  and  gliding  can  be  effectively  stabilized  by 
smooth  wing  motions  driven  by  the  central  pattern  generator  network.  This  work,  published  in  the 
AIAA  Journal  of  Guidance,  Control,  and  Dynamics,  was  awarded  the  Best  Paper  Award  from  the 
AIAA  Infotech  @  Aerospace  Conference  in  2009  [Jl].  Another  contribution  of  this  AFOSR- funded 
research  is  to  investigate  at  the  transient  behaviors  between  flapping  and  gliding  motions  of  birds. 
In  contrast  with  small  insects  or  humming  birds,  birds  do  glide  and  perch  by  controlling  their 
articulated  wings.  New  insights  have  been  derived  into  how  birds  or  bats  without  a  vertical  tail 
can  stabilize  and  control  their  agile  motions  [J2, J3] .  This  part  of  research  is  essential  to  realize 
a  bat-sized  MAV  that  has  to  alternate  between  flapping  and  gliding  for  lower-energetic  costs  and 
perching  maneuvers.  In  contrast  with  other  ongoing  research  projects  in  flapping  flight  of  an  insect¬ 
like  MAV,  flapping  flying  MAVs  whose  size  are  similar  to  that  of  birds  and  bats  are  developed  for 
this  project,  with  an  aim  to  demonstrate  closed-loop  experimentation  of  engineered  flapping  flight. 
A  robotic  bat  with  8  motors  has  been  developed  and  installed  on  a  3-DOF  pendulum  [Bl].  This 
AFOSR  project  successfully  demonstrated  that  complicated  flapping  dynamic  motions  could  be 
controlled  by  using  phase  synchronization  of  CPG  oscillators  by  using  this  robotic  bat  testbed. 

Dynamics  and  Performance  of  Articulated- Wing  Aircraft.  The  dynamics  of  articulated- 
winged  MAVs  without  a  vertical  tail  has  been  studied  in  detail  [J2, J3,B1] .  This  work  shows  that 
how  birds  can  effectively  control  its  dihedral  angles  in  gliding  flight  to  stabilize  their  lateral  motions 
without  a  vertical  tail.  Although  the  dynamics  and  control  of  articulated  wing  aircraft  share  several 
common  features  with  conventional  fixed  wing  aircraft,  the  presence  of  wing  articulation  presents 
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several  unique  benefits  as  well  as  limitations  from  the  perspective  of  performance  and  control.  One 
of  the  objectives  of  this  AFOSR  project  is  to  understand  these  features  using  a  combination  of  the¬ 
oretical  and  numerical  tools.  The  aircraft  concept  envisioned  in  this  project  uses  the  wing  dihedral 
angles  for  longitudinal  and  lateral-directional  control.  Aircraft  with  flexible  articulated  wings  are 
also  investigated.  A  complete  nonlinear  model  of  the  flight  dynamics  incorporating  dynamic  Center 
of  Gravity  (CG)  location  and  the  changing  moment  of  inertia  has  been  derived.  It  is  shown  that 
symmetric  dihedral  configuration,  along  with  a  conventional  horizontal  tail,  can  be  used  to  control 
flight  speed  and  flight  path  angle  independently  of  each  other.  This  characteristic  is  very  useful 
for  initiating  an  efficient  perching  maneuver.  Experimental  results  [B1,J5]  are  presented  which 
help  demonstrate  the  capability  of  wing  dihedral  for  control  and  for  executing  maneuvers  such 
as  slow,  rapid  descent  and  perching.  Both  open-loop  and  closed-loop  experiments  with  the  Cgull 
articulated-wing  MAV  were  performed  to  demonstrate  (a)  the  effectiveness  of  symmetric  dihedral 
for  flight  path  angle  control,  (b)  yaw  control  using  asymmetric  dihedral,  and  (c)  the  elements  of 
perching. 

Dynamics  and  Control  of  Aircraft  with  Flexible,  Articulated  Wings.  It  is  shown  that 
wing  dihedral  angles  alone  can  effectively  regulate  sideslip  during  rapid  turns  and  generate  a  wide 
range  of  equilibrium  turn  rates  while  maintaining  a  constant  flight  speed  and  regulating  sideslip. 
We  compute  the  turning  performance  limitations  that  arise  due  to  the  use  of  wing  dihedral  for  yaw 
control,  and  compare  the  steady  state  performance  of  rigid  and  flexible- winged  aircraft.  We  present 
an  intuitive  but  very  useful  notion,  called  the  effective  dihedral,  which  allows  us  to  extend  some  of 
the  stability  and  performance  results  derived  for  rigid  aircraft  to  flexible  aircraft.  In  the  process,  we 
identify  the  extent  of  flexibility  needed  to  induce  substantial  performance  benefits,  and  conversely 
the  extent  to  which  results  derived  for  rigid  aircraft  apply  to  a  flexible  aircraft.  We  demonstrate, 
interestingly  enough,  that  wing  flexibility  actually  causes  deterioration  in  the  maximum  achievable 
turn  rate  when  the  sideslip  is  regulated. 

PDE  Boundary  Control  of  Flexible  Aircraft.  This  AFOSR  project  also  reports  on  the 
first  PDE  control  result  for  controlling  flexible  wings  of  MAVs  [J4] .  Using  a  simple  order  of  magni¬ 
tude  analysis,  we  derive  conditions  under  which  the  wing  is  structurally  statically  stable,  as  well  as 
conditions  under  which  there  exists  time  scale  separation  between  the  bending  and  twisting  dynam¬ 
ics.  We  show  that  the  time  scale  separation  depends  on  the  geometry  of  the  wing  cross  section,  the 
Poisson’s  ratio  of  the  wing  material,  the  flight  speed  and  the  aspect  ratio  of  the  wing.  We  design 
independent  control  laws  for  bending  and  twisting.  A  key  contribution  of  this  AFOSR  project  is 
the  formulation  of  a  partial  differential  equation  (PDE)  boundary  control  problem  for  wing  de¬ 
formation.  PDE-backstepping  is  used  to  derive  tracking  and  exponentially  stabilizing  boundary 
control  laws  for  wing  twist  which  ensure  that  a  weighted  integral  of  the  wing  twist  (net  lift  or  the 
rolling  moment)  tracks  the  desired  time- varying  reference  input.  We  show  that  a  control  law  which 
only  ensures  tracking  of  a  weighted  integral  improves  the  stability  margin  of  the  twisting  dynamics 
sixteen  fold.  A  tracking  control  law  is  derived  for  the  wing  tip  displacement  which  uses  motion 
planning  and  a  novel  two-stage  perturbation  observer.  This  work  on  PDE-based  control  of  wing 
deformation  allows  for  the  use  of  highly  flexible  wings  on  MAVs. 

Put  together,  this  AFOSR  YIP  project  provides  a  comprehensive  understanding  of  the  flight 
dynamics  of  a  robotic  aircraft  or  MAV  equipped  with  articulated  wings,  and  provides  a  set  of  con- 
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trol  laws  for  performing  agile  maneuvers  and  for  honing  the  benefits  of  using  highly  flexible  wings. 


This  AFOSR  project  has  produced  1  PhD  thesis  and  1  MS  thesis,  while  1  PhD  thesis  is  in 
progress.  Also,  the  PI  has  published  3  journal  papers  and  2  pending  journal  articles.  One  book 
chapter  article  will  be  published  in  the  AIAA  Progress  in  Aeronautics  and  Astronautics  series. 
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1  Introduction 

There  is  a  considerable  interest  in  developing  robotic  aircraft  inspired  by  birds  and  bats  [44,  14,  31] 

and  insects  [24,  22,  21,  107].  In  contrast  with  insects  whose  wings  are  well  modeled  as  simple 
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rigid  wings,  both  wing  flexibility  and  wing  articulation  are  believed  to  play  a  key  role  in  flight 
performance  and  agility  for  bird  and  bat  flight  [90].  The  broader  goal  is  to  learn  and  mimic  avian 
flight  with  the  ultimate  objective  of  designing  unmanned  aerial  vehicles  which  are  autonomous,  agile 
and  capable  of  flying  in  constrained  environments  [57].  Birds  are  natural  role  models  for  designing 
robotic  micro  air  vehicles  (MAVs)  wherein  the  aforementioned  attributes  can  be  engineered.  MAVs 
typically  fly  in  a  low  Reynolds  number  range  of  103  —  105  [54]  which  coincides  with  that  of  the 
birds.  Therefore,  it  is  worth  investigating  the  mechanics  of  avian  flight  and  making  an  attempt 
to  reverse-engineer  them.  Conversely,  a  study  of  the  flight  mechanics  of  MAVs  can  shed  light  on 
several  aspects  of  bird  flight. 

Birds  lack  a  vertical  tail.  Instead,  the  wing  dihedral  and  incidence  angles  are  controlled  actively 
as  birds  execute  agile  and  even  spatially  constrained  maneuvers.  Birds  are  known  to  have  a  very 
complex  wing  structure,  with  a  wide  range  of  “actuators.”  The  wings  can  flap,  twist  and  change 
the  sweep  angle  on  demand.  They  have  a  wide  range  of  feathers  which  serve  as  flaps  and  spoilers. 
The  hair  on  bird  wings  can  sense  local  flow  conditions,  and  feedback  from  these  sensors  is  sent  to 
the  feathers  which  are,  either  passively  or  actively,  oriented  to  optimize  the  flow  conditions  on  the 
wing  for  stability  and  maneuverability. 

Over-actuation  to  the  extent  seen  in  birds  is  difficult  to  engineer  in  aircraft  because  of  well  un¬ 
derstood  considerations  such  as  weight  penalty,  actuator  limitations,  sensor  design,  etc.  Therefore, 
it  is  necessary  to  abstract  out  the  underlying  phenomena  and  understand  their  implications  for 
stability  and  control.  That  way,  they  can  be  engineered  onto  actual  aircraft.  This  is  the  motivation 
behind  this  AFOSR  project. 


Deformable 

wings 


(a)  Robotic  aircraft  concept  (b)  Robotic  bat  at  UIUC 

Figure  1:  Figure  (a)  shows  a  robotic  aircraft  with  flexible  wings  controlled  at  the  wing  root  (by 
servo  actuators)  and  the  tip  (by  flaps).  Figure  (b)  shows  a  robotic  bat  testbed  where  the  control 
laws  proposed  in  this  project  can  be  tested  [44,  14], 

This  AFOSR  project  contributes  to  the  broader  problem  of  developing  a  flapping  MAV  capable 
of  agile  flight  in  constrained  environments,  such  as  the  conceptual  MAV  in  Fig.  1(a).  Chung  and 
Dorothy  [14]  studied  a  neurobiologically-inspired  controller  for  flapping  flight,  and  demonstrated  it 
on  a  robotic  testbed  shown  in  Fig.  1(b).  Their  controller  could  switch  in  a  stable  and  smooth  fashion 
between  flapping  and  gliding  flight.  Gliding  is  essential  during  soaring,  landing  and  perching  (see 
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Fig.  2  for  a  schematic  diagram  of  a  perching  maneuver).  This  project  aims  to  provide  the  analytical 
foundations  for  using  wing  articulation  to  perform  agile  maneuvers  while  gliding  and  landing. 


Figure  2:  Schematic  of  a  perching  maneuver  executed  by  a  tailless  articulated  wing  aircraft. 

Like  any  problem  in  flight  mechanics  and  control,  this  project  approaches  the  chosen  problem 
along  three  directions:  performance,  stability  and  control.  Each  facet  has  been  addressed  in  the 
report: 

•  Performance:  we  focus  on  steady  state  performance  metrics  such  as  the  turn  rate,  gliding 
angle  and  speed.  Turn  rate  is  of  particular  interest  because  this  performance  metric  would  be 
most  affected  by  the  absence  of  a  vertical  tail.  We  also  analyse  the  influence  of  wing  flexibility 
on  the  performance  metrics. 

•  Stability:  we  use  bifurcation  analysis  to  study  the  nature  of  stability  of  the  turn  and  wings 
level  equilibria.  We  derive  literal  approximations  to  some  lateral  stability  indicators. 

•  Control:  we  demonstrate  control  and  control-related  issues  through  a  combination  of  theoret¬ 
ical  analysis  and  experiments  on  an  actual  MAV.  We  also  design  PDE-based  control  schemes 
to  control  the  deformation  of  a  flexible  wing.  The  control  design  process  makes  use  of  several 
features  of  MAV  geometry  and  flight  regime. 

One  point  that  is  particularly  worth  highlighting  is  the  restriction  of  our  analysis  to  MAVs. 
This  restriction  comes  about  due  to  the  fact  that  wing  articulation  proposed  in  this  project  would 
be  impractical  for  larger  aircraft.  The  low  flight  speed  and  small  size  of  MAVs  lead  to  several 
interesting  consequences.  The  low  flight  speeds  imply  that  the  local  angle  of  attack  distribution  on 
the  wing  is  a  nonlinear  function  of  the  angular  velocity  of  the  aircraft.  The  nonlinearity  manifests 
itself  in,  for  example,  the  peculiar  relation  between  the  turn  rate  and  the  commanded  dihedral  as 
described  later  in  the  report.  The  low  flight  speed  and  the  Young’s  modulus  together  lead  to  time 
scale  separation  between  the  twisting  and  bending  dynamics  of  a  flexible  wing,  which  simplifies  the 
control  design  process  significantly  in  practical  situations. 

The  problems  that  we  have  considered  in  the  AFOSR  project  have  several  novel  elements. 
Before  stating  the  contributions  of  the  project,  we  briefly  discuss  the  state  of  the  art. 
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2  Literature  Review 


Fixed  and  flapping  wing  MAVs  have  been  extensively  studied  in  the  literature.  The  reader  is 
referred  to  an  excellent  compendium  of  papers  [54]  which  showcases  some  of  the  work  done  in 
this  area  until  circa  2000.  The  idea  of  using  cant-angle  winglets  for  control  was  investigated 
experimentally  by  Bourdin,  Gatto  and  Friswell  [8,  28]  for  larger  aircraft  operating  at  high  Reynolds 
numbers.  However,  their  analysis  was  restricted  to  calculating  the  aerodynamic  moments  arising 
from  the  use  of  cant-angle  winglets.  Wickenheiser  and  Garcia  [105,  106]  studied  the  dynamics  of 
morphing  aircraft  and  demonstrated  perching  using,  among  other  forms  of  articulation,  variable 
wing  incidence.  Reich  et  al.  [76]  experimentally  studied  the  aerodynamic  performance  of  a  wing  of 
variable  incidence  for  perching.  Crowther  [18]  showed  that  a  perched  landing  can  be  achieved  using 
a  simple  pitch-up  maneuver,  and  a  similar  conclusion  was  obtained  after  optimization  by  Cory  and 
Tedrake  [17]. 

Stenfelt  and  Ringertz  [92,  93]  studied  the  stability  and  control  of  tailless  aircraft  equipped  with  a 
split  flap  mechanism.  Shtessel,  Buffington  and  Banda  [86]  designed  a  sliding  mode-based  controller 
for  tailless  fighter  aircraft,  while  Patel  and  co-authors  [66]  designed  a  robust  adaptive  controller 
for  tailless  aircraft  in  the  presence  of  actuator  failures.  Recently,  Obradovic  and  Subbarao  [56] 
investigated  the  power  requirement  for  morphing,  and  used  it  as  a  basis  to  compare  wing  morphing 
and  the  traditional  aileron-based  control  in  different  flight  regimes.  The  lateral  stability  and  control 
of  birds,  and  in  particular,  the  role  of  wing  dihedral,  have  been  studied  extensively  by  Sachs  and 
co-authors  [81,  79,  80].  Sachs  has  demonstrated  that  for  air  vehicles  whose  size  and  speed  (and 
hence,  the  Reynolds  number)  are  similar  to  those  of  birds,  wings  are  sufficient  to  provide  lateral 
stability  thereby  reducing,  if  not  eliminating  altogether,  the  need  for  a  vertical  tail.  Tran  and  Lind 
[99]  numerically  computed  the  stability  of  an  aircraft  equipped  with  variable  symmetric  dihedral 
and  incidence.  Their  wing  model  consisted  of  an  assemblage  of  rigid  segments. 

A  variety  of  aircraft  models  incorporating  wing  and  fuselage  flexibility  have  been  proposed  in  the 
literature,  although  most  of  these  models  do  not  consider  wing  articulation.  Waszak  and  Schmidt 
[104]  derived  a  complete  nonlinear  model  of  an  aircraft  with  flexible  wings.  Their  aerodynamic 
model,  however,  assumed  a  steady  flow,  and  their  frame  of  reference  consisted  of  the  so-called 
mean  axes  which  are  hard  to  locate  in  a  practical  situation.  Tuzcu  and  Meirovitch  [52]  extended 
their  model  in  several  ways:  they  used  a  more  intuitive  reference  frame  (the  conventional  body 
axes)  and  a  more  accurate  Theodorsen’s  unsteady  aerodynamics  theory  for  computing  the  forces 
and  moments  [97].  Recently,  Nguyen  and  Tuzcu  [55]  presented  a  dynamic  model  for  a  fully  flexible 
aircraft.  These  papers  worked  with  a  small-strain,  small-displacement  beam  theory.  In  contrast, 
Patil  and  co-authors  [67,  71]  derived  a  geometrically  exact  (large  displacement)  small-strain  nonlin¬ 
ear  beam  model,  and  used  it  to  study  the  dynamics  and  stability  of  flying  wings.  Shearer  and  Cesnik 
[85]  and  Su  and  Cesnik  [96]  used  nonlinear  flight  dynamic  and  structural  models  to  investigate  the 
effects  of  structural  nonlinearities  on  the  dynamic  stability  of  aircraft  characterized  by  large  aspect 
ratio  wings  and  blended  wing-body  configurations,  respectively.  Baghdadi  and  co-authors  [3]  used 
bifurcation  analysis  to  study  the  performance  and  stability  of  a  flexible  aircraft  model  based  on 
Ref.  [104]  and  concluded  that  flexibility  must  be  accounted  for  carefully  during  the  control  design 
process.  They  also  demonstrated  that  a  control  law  designed  assuming  a  rigid  configuration  could 
trigger  instabilities  in  an  otherwise  identical  aircraft  with  flexible  wings.  Rodden  [77]  derived  ana¬ 
lytical  expressions,  backed  by  experimental  approximations,  for  increments  in  the  rolling  moment 
derivatives  arising  from  aeroelastic  effects. 

It  has  been  demonstrated  in  the  literature  that  aeroelastic  instabilities  such  as  wing  divergence 
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and  flutter  can  be  mitigated  using  flap-based  effectors  [5,  74]  or  passive  energy  sinks  (for  flutter) 
[36].  There  is  a  substantial  amount  of  literature  on  boundary  control  theory  of  PDEs  (see  Refs.  [78, 
43,  42,  12,  13,  47]  for  material  pertinent  to  this  report  and  the  references  cited  therein).  There  are 
two  sets  of  methods  for  boundary  control  of  PDEs.  The  first  set  consists  of  methods  that  seek  to 
convert  the  PDEs  into  ordinary  differential  equations  using  approximation  methods  such  as  those 
of  Galerkin  or  Rayleigh- Ritz  [13,  34],  or  using  operator  theory  [19,  9].  The  second  set  consists  of 
methods  that  keep  the  PDEs  intact,  and  use  a  “model-following”  approach  as  described  in  a  recent 
book  by  Krstic  and  Smyshlyaev  [42] . 

If  a  PDE  is  approximated  by  Galerkin’s  method,  or  converted  into  an  ODE  form  using  operator 
theory,  the  problem  of  achieving  an  integral  objective  reduces  to  a  standard  output  control  problem. 
Whereas  solutions  to  output  control  problems  in  an  ODE  setting  are  abundantly  known,  an  ODE- 
based  approximation  to  PDEs  usually  leads  to  systems  having  large  orders.  The  ODE-based  control 
design  process  usually  becomes  tedious,  the  control  laws  become  non-intuitive,  and  the  tracking 
errors  could  be  large  for  a  class  of  inputs  when  a  poor  choice  is  made  for  the  basis  functions 
which  capture  the  time-varying  boundary  conditions.  Ref.  [12]  also  points  out  that  a  finite-state 
approximation  may  wrongly  render  fundamental  system  theoretic  properties  like  controllability 
and  observability  to  be  functions  of  the  approximation.  Stability  analysis  based  on  finite  state 
approximation  is  vulnerable  to  spillover  instabilities  which  arise  due  to  inadequate  accounting  of 
the  residual  modes  [4,  51]. 

On  the  other  hand,  keeping  PDEs  intact  makes  the  control  law  design  more  intuitive.  It  has 
been  used  in  the  past  for  maneuvering  robotic  arms  [70,  42],  controlling  the  Navier-Stokes  model 
[102]  and  suppressing  vibrations  in  a  flexible  beam  [32],  A  gain-scheduling  based  approach  for 
nonlinear  PDEs  has  been  presented  in  Ref.  [88],  while  Krstic  and  Smyshlyaev  [43]  derived  an 
adaptive  controller  for  parabolic  PDEs. 

3  Main  Contributions 

3.1  Performance  and  Stability 

As  we  stated  earlier,  this  AFOSR  project  is  dedicated  to  MAVs  whose  speed  and  size  present  several 
distinct  characteristics  from  the  point  of  view  of  control  and  stability.  The  small  size  of  MAV  wings 
makes  wing  articulation  practically  feasible.  Unlike  conventional  fixed  wing  aircraft,  an  articulated 
wing  aircraft  changes  its  configuration  routinely  and,  therefore,  stability  is  closely  tied  to  the  nature 
of  the  maneuver  being  executed.  We  use  bifurcation  analysis  [30,  48,  65]  to  explore  the  dynamics 
of  tailless  aircraft  equipped  with  articulated  wings.  Bifurcation  analysis  not  only  measures  the 
stability  characteristics  of  the  aircraft,  but  also  helps  predict  the  performance  limitations  that  arise 
because  of  the  use  of  asymmetric  dihedral.  We  study  a  rigid  aircraft  model  to  gain  a  foothold,  and 
then  perform  a  similar  analysis  on  an  aircraft  with  flexible  wings  for  comparison. 

The  report  includes  detailed  theoretical  and  linear  computational  analyses  of  the  lateral  dy¬ 
namics.  Longitudinal  dynamics  are  not  affected  by  the  absence  of  a  vertical  stabilizer.  Analytical 
expressions  for  lateral-directional  aerodynamic  force  and  moment  derivatives  offer  valuable  insight 
into  the  maneuver-dependence  of  stability  and  help  identify  the  source  of  lateral-directional  insta¬ 
bility,  which  is  subsequently  verified  computationally.  The  analytical  expressions  also  help  identify 
potentially  dangerous  situations  where  the  control  effectiveness  of  the  dihedral  may  switch  sign  in 
the  midst  of  certain  maneuvers. 


It  is  known  that  the  wing  dihedral  angle  can  be  varied  to  perform  slow,  steep  descents  [81,  98]. 
We  compute  the  gliding  flight  equilibria  numerically  along  with  their  stability  to  identify  bounds 
on  longitudinal  performance.  A  knowledge  of  the  longitudinal  trims  can  help  formulate  landing 
and  perching  strategies  in  spatially-constrained  environments  without  resorting  to  maneuvers  like 
spin  [75]  and  aid  the  design  of  control  laws  for  perching  [10]. 

Next,  we  compare  select  performance  metrics  and  stability  of  the  rigid  MAV  with  those  of  a 
similar  MAV  equipped  with  flexible  wings.  The  purpose  of  this  analysis  is  two-fold.  First,  it  helps 
to  identify  the  benefits  and  limitations  of  using  wing  flexibility.  Second,  and  less  obviously,  it  helps 
to  identify  the  extent  to  which  a  rigid  MAV  model  can  accurately  capture  the  dynamics  of  a  flexible 
wing  MAV.  The  flight  dynamics  will  be  rendered  unstable  if  the  wing  is  structurally  unstable.  On 
the  other  hand,  if  the  structural  stability  of  the  wing  can  be  guaranteed,  then  the  performance 
and  stability  of  the  motion  can  be  computed  reliably  by  considering  “macroscopic”  parameters  like 
the  resultant  forces  which  depend  on  the  shape,  rather  than  stability,  of  the  wing.  Therefore,  an 
analysis  like  the  one  in  this  project  would  depend  largely  on  the  wing  geometry  (which  is  usually 
well  known  a  priori )  rather  than  a  precise  knowledge  of  the  elastic  parameters. 

The  specific  questions  answered  for  a  flexible-winged  MAV  include: 

1.  For  the  wing  of  a  typical  bird-sized  MAV,  what  value  of  Young’s  modulus  (E)  should  the  wing 
have  in  order  for  the  MAV  to  offer  a  significant  performance  improvement  over  a  rigid  wing 
MAV?  Equivalently,  until  what  point  is  the  open  loop  analysis  of  a  rigid  aircraft  relevant  to 
a  flexible-winged  aircraft?  The  notion  of  effective  dihedral  is  introduced  in  a  bid  to  answer 
these  questions. 

2.  A  stiff  wing  may  be  required  for  certain  maneuvers.  Axial  tension  in  the  wing  is  an  intuitive 
stiffener.  How  effective  and  useful  is  it?  We  answer  this  question  in  the  negative  -  it  is 
effective,  but  of  limited  utility. 

3.  How  is  the  stability  of  the  motion  altered  in  the  presence  of  flexible  wings?  The  wings  are 
assumed  to  be  quasi-statically  deformed  and,  therefore,  the  structural  dynamics  of  the  wings 
have  no  bearing  on  the  conclusion.  In  other  words,  the  wing  is  assumed  to  be  structurally 
stable  and  its  dynamics  sufficiently  faster  than  the  aircraft. 

4.  Is  there  a  measurable  improvement  in  the  steady  state  turning  performance?  Steady  state 
turn  rate  is  the  only  agility  metric  which  is  based  entirely  on  a  steady  maneuver  [61].  It  is 
also  an  important  benchmark  to  evaluate  the  efficacy  of  a  yaw  control  mechanism. 

3.2  Control 

The  ideas  discussed  above  were  subjected  to  experimental  validation  involving  open  loop  as  well  as 
closed  loop  flight  tests.  The  purpose  of  the  experiments  is  to  demonstrate  (a)  the  effectiveness  of 
symmetric  dihedral  for  flight  path  angle  control,  (b)  yaw  control  using  asymmetric  dihedral,  and 
(c)  some  elements  of  a  perching  maneuver  [26]. 

One  of  the  key  contributions  of  this  AFOSR  project  is  a  boundary  control  problem  for  wing  twist 
which  could  be  extended  to  a  wider  class  of  hyperbolic  partial  differential  equations  (PDEs).  We 
also  consider  motion  planning  for  the  bending  dynamics  of  an  Euler-Bernoulli  beam  using  a  novel 
estimator-based  scheme.  The  problem  of  controlling  a  flexible  wing  presents  unique  problems  as 
well  as  helpful  features.  A  flexible  wing  represents  a  coupled  twist-bending  problem  in  the  simplest 
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case,  and  a  coupled  twist-bending-rigid-motion  problem  in  the  most  general  case.  Either  way, 
problems  with  coupled  PDEs  are  difficult  to  solve  because  the  only  unified  approach  to  assembling 
a  system  of  diverse  PDEs  requires  that  they  be  converted  into  an  equivalent  or  approximate  ODE 
system.  This  is  a  cumbersome  route. 

We  show  that  the  twisting  dynamics  of  the  MAV  wing  are  faster  than  the  bending  dynamics, 
and  the  time  scale  separation  depends  on  the  aspect  ratio  of  the  wing.  We  design  independent 
tracking  controllers  for  twist  and  bending,  which  can  be  applied  effectively  to  control  the  coupled 
twist-bending  dynamics.  We  also  identify  the  extent  of  flexibility  required  to  trigger  a  structural 
instability  in  an  MAV  at  routine  flight  speeds. 

We  show  that  wing  twist  can  be  stabilized  rapidly  using  a  root-based  actuator  driven  by  a  PDE 
backstopping  controller.  The  procedure  is  called  backstepping  because  it  involves  a  Volterra  oper¬ 
ator  with  a  lower  triangular  structure  similar  to  backstepping  transforms  for  ordinary  differential 
equations  [41,  42].  It  is  a  continuum  analogue  of  the  backstepping  transformations  in  ODEs  and 
allows  the  controller  acting  at  the  boundary  to  compensate  for  the  undesired  (unmatched)  dynam¬ 
ics.  We  also  design  a  tracking  control  (which  is  considerably  harder  when  the  actuator  is  located 
only  at  the  boundary)  based  on  motion-planning  which  can  be  added  on  top  of  the  backstepping 
controller. 

For  a  class  of  input-output  combinations,  we  show  that  the  twist  dynamics  have  a  finite  relative 
degree  which  permits  a  more  traditional  approach  to  control  design.  We  show  that  the  tracking 
controller  renders  the  cantilever  wing  into  the  form  of  a  clamped-clamped  beam,  which  improves  the 
stability  margin  of  the  twist  dynamics  by  a  factor  of  sixteen.  An  adaptive  controller  for  a  limited 
class  of  parametric  uncertainties  is  also  derived  for  tip-based  actuators  in  general,  and  root-based 
actuators  when  the  output  is  the  rolling  moment. 

There  are  several  controller  designs  in  the  literature  for  stabilizing  and  controlling  beam  bending 
(see,  for  example,  [32,  88,  101,  27,  40]).  We  design  a  perturbation  observer-based  controller  to 
facilitate  a  motion  planning-based  tracking  controller  for  bending.  The  output  of  interest  is  the 
displacement  of  the  wing  tip.  As  the  name  suggests,  the  perturbation  observer  designed  here 
uses  adaptation  to  estimate  the  external  forces  acting  on  the  system.  The  observer  is  split  into 
a  “particular”  and  a  “homogeneous”  component  (the  notions  are  made  more  precise  later).  Since 
the  homogeneous  component  is  stable  and  not  driven  directly  by  external  feedback,  it  is  simpler  to 
design  a  control  law  for  it.  The  same  control  signal  is  sent  to  the  actual  system,  whose  states  then 
converge  exponentially  rapidly  to  a  bounded  envelope  around  the  observer  states. 

4  Biologically  Inspired  CPG-Based  Control 

The  flapping  motion  in  bats  combines  many  different  joints  and  muscles  to  create  rhythmic  motions 
which  produce  the  aerodynamic  forces  required  to  sustain  flight.  With  the  current  state  of  the 
actuator /structure  technology,  we  are  unable  to  mimic  the  strength  and  flexibility  of  the  distributed 
structure  of  bones,  joints  and  ligaments.  More  specifically,  we  are  unable  to  produce  muscle- 
strength  actuation  of  the  entire  system  using  mechanical  actuators.  For  this  purpose,  we  have 
settled  on  starting  our  investigation  with  three  main  modes  of  motion:  flapping,  twisting,  and 
lead-lag  (sometimes  referred  to  as  plunging,  pitching,  and  sweeping).  These  three  motions  (shown 
in  Figure  3)  have  the  strongest  connection  to  traditional  flight  models  and  were  the  the  first  to  be 
discussed  in  observations  of  animal  flight  [2].  Azuma  also  notices  that  the  phase  differences  between 
the  motions  is  extremely  important  [2].  Next  we  consider  how  to  use  CPGs  for  precise  control  of 


10 


(a)  RoboBat  Flapping  Motion, 

<t> 


(b)  RoboBat  Pitching  Motion, 

e 


(c)  RoboBat  Lead-Lag  Motion, 


Figure  3:  Main  wing  motions  used  by  RoboBat 


phase  differences  between  the  three  motions. 

4.1  Central  Pattern  Generators 

Many  creatures  produce  their  motion  by  synchronizing  periodic  motions  of  limbs,  such  as  running, 
swimming  or  flapping.  They  do  this  by  coupling  biological  oscillators  and  synchronizing  their  out¬ 
puts.  Biological  oscillators  rely  on  short  timescale  (ms)  neuron  dynamics  including  spike-bursting, 
spike  frequency  adaptation,  and  post-inhibitory  rebound.  Herrero-Carron,  et.  al.  [33]  designed  a 
control  law  for  modular  robots  by  approximating  short  timescale  neuron  dynamics.  Because  there 
is  such  a  short  timescale  required  for  integration,  the  neuron  dynamics  were  integrated  offline.  We 
are  unlikely  to  be  able  to  perform  such  strict  mimicry  in  an  online  controller  as  we  add  additional 
neurons  for  feedback,  active  control  of  phase  differences,  or  gait  transitions. 

In  order  to  make  online  control  more  feasible,  we  can  emulate  these  biological  oscillators  by  using 
coupled  nonlinear  limit  cycle  oscillators.  A  limit  cycle  oscillator  converges  to  a  stable  trajectory 
which  is  called  the  limit  cycle.  Because  of  this  convergence  the  oscillator  will  quickly  forget  sporadic 
disturbances  and  converge  back  to  the  stable  limit  cycle.  If  the  oscillator  itself  is  a  smooth  vector 
field,  we  can  smoothly  transition  between  desired  trajectories  without  abrupt  changes  being  required 
in  the  motor  output. 

This  type  of  control  design  allows  the  complicated  synchronization  and  trajectory  computations 
to  be  performed  according  to  simple  rules  that  provide  the  desired  oscillatory  behavior.  If  we  can 
then  influence  body  motion  by  simply  tuning  top-level  inputs  into  the  CPG  network,  we  can  achieve 
the  dimension  reduction  that  allows  the  control  problem  to  be  computationally  feasible.  Vertebrate 
brains  send  top-level  chemical  signals  which  modulate,  start,  and  stop  CPG  behavior,  but  do  not 
micromanage  the  joint  trajectories  [91]. 

Following  our  previous  work  [14],  we  use  the  following  limit-cycle  model  called  the  Hopf  oscil¬ 
lator,  named  after  the  supercritical  Hopf  bifurcation  model  with  a  =  1: 

and  (1) 


where  the  A  >  0  denotes  the  convergence  rate  to  the  symmetric  limit  circle  of  the  radius  p  >  0 


x  =  f(x;  p ;  cr)  +  u(t) 
-A 


with  x  =  (u  —  a,  v)J 


f(x;p;cr)  = 


(u-a)2+v2 


-A 


-u(t) 

(u-a)2-\-v2 
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(a)  Configuration  A 


(b)  Symmetric  Configuration  A 


(c)  Configuration  B 


Figure  4:  Example  graph  configurations  of  the  coupled  Hopf  oscillators  on  balanced  graphs. 


and  u (t)  is  an  external  or  coupling  input.  The  bifurcation  parameter  a  can  change  to  a  negative 
number  (e.g.  —1  such  that  — b  1^).  This  would  change  the  stable  limit  cycle  dynamics 

to  the  dynamics  with  a  globally  stable  equilibrium  point  at  the  bias  “a”  [95].  Such  a  change  can 
be  used  to  turn  the  flapping  oscillatory  motion  to  the  gliding  mode,  as  was  seen  in  [14].  This 
could  be  used  to  mimic  many  animal  flight  strategies.  Later  sections  of  this  paper  discuss  terminal 
perching  control,  assumed  to  be  in  a  gliding  phase.  Often,  bats  and  birds  will  transition  to  gliding 
in  order  to  shed  energy  well  before  an  actual  perching  maneuver  is  performed.  In  [14],  it  was  shown 
that  coupled  networks  of  Hopf  oscillators  on  balanced  graphs  exhibit  smooth  exponentially  stable 
behavior  in  both  oscillatory  mode  and  fixed  point  mode. 

The  possibly  time- varying  parameter  uj(t)  >  0  determines  the  oscillation  frequency  of  the 
limit  cycle.  A  time-varying  a{t)  sets  the  bias  to  the  limit  cycle  such  that  it  converges  to  u(t)  = 
p  cos  (c ot  +  (5)  +  a  and  v(t)  =  psin  (c ot  +  <5)  on  a  circle.  The  output  variable  to  generate  the  desired 
oscillatory  motion  of  each  joint  is  the  first  state  u  from  the  Hopf  oscillator  model  in  Eq.  (1). 

Phase  synchronization  means  an  exact  match  of  the  scaled  amplitude  or  the  frequency  in  this 
paper.  Hence,  phase  synchronization  permits  different  actuators  to  oscillate  at  the  same  frequency 
but  with  a  prescribed  phase  lag.  In  essence,  each  CPG  dynamic  model  in  Eq.  (1)  is  responsible 
for  generating  the  limiting  oscillatory  behavior  of  a  corresponding  joint,  and  the  diffusive  coupling 
among  CPGs  reinforces  phase  synchronization.  For  example,  the  flapping  angle  has  roughly  a  90- 
degree  phase  difference  with  the  pitching  joint  to  maintain  the  positive  angle  of  attack  (see  the 
actual  data  from  birds  in  [2]).  The  oscillators  are  connected  through  diffusive  couplings,  and  the 
i-th  Hopf  oscillator  can  be  rewritten  with  a  diffusive  coupling  with  the  phase-rotated  neighbor. 


±i  =  f(x;;  Pi)  -  k  x  ^ 


E 

jeMi 


Xj  -  —  R(Ajj)xj 

Pi 


(2) 


where  the  Hopf  oscillator  dynamics  f  (x?; ;  p{)  with  a  =  1  is  defined  in  Eq.  (1),  J\fi  denotes  the  set 
that  contains  only  the  local  neighbors  of  the  i-th  Hopf  oscillator,  and  m,  is  the  number  of  the 
neighbors.  The  2x2  matrix  R(Ay)  is  a  2-D  rotational  transformation  of  the  phase  difference  Ay 
between  the  i-th  and  j-th  oscillators.  The  positive  (or  negative)  Ay  indicates  how  much  phase  the 
i-th  member  leads  (or  lags)  from  the  j-th  member  and  Ay  =  —A The  positive  scalar  k  denotes 
the  coupling  gain. 

Numerous  configurations  are  possible  as  long  as  they  are  on  balanced  graphs  [16]  and  we  can 
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choose  either  a  bidirectional  or  a  uni-directional  coupling  between  the  oscillators.  The  numbers 
next  to  the  arrows  in  Fig.  4  indicate  the  phase  shift  A y,  hence  AtJ  >  0  indicates  how  much  phase 
the  i-tlr  member  leads.  Figure  4(a)  shows  the  choice  of  coupling  used  in  this  paper,  where  <f)w ,  9W, 
and  ipw  are  flapping  angle,  twist  angle,  and  lead-lag  angle,  respectively.  Subscripts  L  and  R  refer 
to  left  and  right  wings.  The  graph  is  balanced.  Further,  all  the  phase  shifts  (A.,j)  along  one  cycle 
add  up  to  a  modulo  of  2ir.  The  proof  of  stability  of  (2)  is  given  in  [14],  As  mentioned  in  [25],  this 
type  of  oscillator  generalizes  other  types  of  waveform  control  such  as  the  split-cycle  [23]. 

4.2  Control  Design  from  Physical  Intuition  and  Biological  Observation 

Our  previous  work  [14]  provided  a  simple  intuition  for  phase  difference  control  of  longitudinal 
motion  in  flapping  flight.  With  a  zero  bias  lead-lag  and  a  center  of  gravity  coinciding  with  the 
stroke  plane,  a  phase  difference  of  270  deg  between  the  flapping  CPG  and  the  lead-lag  CPG  gives 
Azuma’s  [2]  elliptical  model  of  flapping:  negative  lead-lag  on  downstroke,  positive  lead-lag  on 
upstroke.  The  simplest  analysis  combines  a  maximum  force  with  the  most-negative  lead-lag  at  the 
middle  of  the  downstroke  to  predict  a  pitch-down  moment  on  the  body.  Alternatively,  if  we  see 
the  phase  difference  to  90  deg,  we  see  the  maximum  force  coinciding  with  the  maximum  positive 
lead-lag  at  the  middle  of  the  downstroke,  predicting  a  pitch-up  moment.  This  intuition  has  been 
confirmed  in  numerical  simulations  and  open  loop  experiments  on  the  RoboBat[14,  15,  44], 


5  Equations  of  Motion  for  an  Articulated  Wing  Aircraft 

The  dynamic  model  derived  in  this  chapter  is  general  enough  that  it  can  be  applied  to  a  wider  class 
of  problems  such  as  flapping  and  a  complete  aeroelastic  analysis  of  aircraft.  This  chapter  consists 
of  two  sections.  The  equations  of  motion  for  a  flexible  aircraft  are  derived  in  the  first  section,  and 
specialized  to  a  rigid  aircraft  in  the  second  section. 


5.1  Notation 


Capital  letters  are  reserved  for  forces,  matrices,  and  for  denoting  coordinate  frames.  Small  letters 
are  used  for  scalars  when  not  in  bold,  and  for  vectors  when  used  with  bold  font.  Given  a  vector 
x  e  R3,  S'(x)  denotes  the  cross  product  operator,  i.e.;  for  any  two  vectors  x,  y  £  R3,  S(x)y  =  xxy. 
Similarly,  52(x)y  =  5(x)(S'(x)y)  =  x  x  (x  x  y).  Given  a  variable  p(t,y),  its  time  derivative  is 


denoted  by  p(t,  y)  = 


dp(t,  y) 


dt 


Its  spatial  derivative  is  denoted  by  p'(t,y )  = 


when  p(t,y)  =p(t),  p(t)  = 


dp(t ) 
dt 


dp(t,  y) 
dy 


Note  that 


5.2  Coordinate  Frames  of  Reference 

Given  frames  F  and  G,  the  matrix  Tpc  is  a  rotation  matrix  which  transforms  the  components  of 
a  vector  from  the  G  frame  to  F.  The  body  frame,  denoted  by  B,  is  attached  to  the  body  with  the 
x  —  z  plane  coincident  with  the  aircraft  plane  of  symmetry  when  the  wings  are  undeflected.  The  x 
axis  points  towards  the  aircraft  nose.  The  z  axis  points  downwards,  and  the  y  is  defined  to  create 
a  right  handed  coordinate  system.  The  coordinate  frames  have  been  illustrated  in  Fig  5,  together 
with  the  dimensions  of  the  aircraft  considered  for  numerical  analysis. 
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26.1  cm 


Figure  5:  A  schematic  of  the  aircraft  showing  the  dimensions  and  the  coordinate  systems  used  to 
model  the  aircraft. 


Consider  the  frame  R  based  at  the  right  wing  root.  Its  origin  coincides  with  that  of  the  B 
frame,  which  is  akin  to  neglecting  the  fuselage  width.  This  assumption  does  not  alter  the  rotation 
matrices  in  any  way.  The  frame  R  is  related  to  the  B  frame  via  three  rotations  at  the  wing  root: 
a  sweep  rotation  Pr  about  the  2-axis,  followed  by  a  dihedral  rotation  Sr  about  the  —x-axis,  and  a 
rotation  Or  about  the  y  axis.  The  y  axis  points  along  the  wing  elastic  axis.  Thus, 


cos  Pr 

-  sin  P R 

0  ' 

'  1 

0 

0 

COS  Or 

0 

sin  Or 

T  RR  = 

sin  pR 

cos  Pr 

0 

0 

COS  Sr 

sin  SR 

0 

1 

0 

0 

0 

1 

0 

-  sin  Sr 

COS  Sr 

_  —  sin  Or 

0 

COS  Or 

A  similar  matrix  can  be  defined  for  the  left  wing.  The  matrix  T rr  is  introduced  here  in  the  most 
general  form,  i.e. ,  no  rotation  is  ignored,  which  makes  it  applicable  to  flapping  flight  dynamics. 
From  hereon,  we  assume  that  Pr  =  Pl  =  0. 

The  frame  S  =  S(y)  is  the  frame  located  at  a  spanwise  wing  station  with  origin  at  the  elastic 
center,  and  y  axis  pointing  along  the  elastic  axis.  The  frame  S  is  related  to  R  via  two  rotations: 
a  rotation  about  the  x  axis  through  the  strain  £'(?/),  and  a  rotation  (twist)  0s(y)  about  the  y  axis. 
Thus, 


Trs 


cos  0s(y)  0  sin  0s(y) 

sm(£'(y))smOs(y)  cos(£'(y))  -  sin(£'(y))  cos  0s{y) 

-cos(f'(y))sin 0s{y)  sin (£'(y))  cos(£'(y))  cos  0s(y) 


(4) 


In  the  interests  of  analytical  tractability,  for  the  purposes  of  computing  the  velocities  and  accel¬ 
eration  terms,  it  will  be  assumed  that  T^g  =  T rr,  i.e.,  the  deformations  are  small  enough  that 
they  do  not  alter  the  coordinate  transformations.  However,  in  Sec.  5.7,  this  assumption  is  relaxed 
for  computing  the  aerodynamic  forces  and  moments.  This  is  the  primary  source  of  the  difference 
in  the  forces  and  moments  produced  by  rigid  and  flexible  wings. 
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5.3  Calculating  the  Velocity  at  a  Spanwise  Station 

The  angular  velocity  of  the  right  wing,  uR  with  components  in  the  body  frame,  is  given  by 

"  -  cos  (3R  -  cos  SR  sin  fjR  0  1  I"  dR 

cjr  =  -  sin  Br  cos  5r  cos  /3r  0  0R  (5) 

0  sin  5r  1  J  [  $r  _ 

It  is  calculated  using  a  3-1-2  Euler  angle  sequence  which  is  also  used  to  calculate  TRr-  The  same 

sequence  can  be  used  to  model  a  flapping  wing,  in  which  case,  the  amplitude  and  phase  of  the 
motion  corresponding  to  each  degree  of  freedom  needs  to  be  prescribed.  In  contrast,  most  flapping 
wing  models  prefer  to  identify  a  stroke  plane  in  which  the  flapping  motion  is  constrained,  and  which 
also  contains  the  twist  axis  (see  Refs.  [14,  7,  87]  and  the  references  cited  therein). 

Let  y  =  [0  y  £]T  denote  the  coordinates  of  a  spanwise  station  on  the  wing  along  the  twisting 
axis.  Then  the  local  wind  velocity,  with  components  in  the  local  station  frame,  is  given  by 

V  =  Tssug  +  u/  +  TsbS(u>b  +  Wi?)TaRy  (6) 

A  similar  expression  can  be  determined  for  the  angular  velocity  of  the  left  wing  at  the  root,  uR, 
and  the  local  velocity  at  a  spanwise  station  on  the  left  wing. 

The  aerodynamic  center  is  assumed  to  be  located  at  the  quarter  chord  point.  The  velocity  at 
the  3/4-chord  point  of  a  spanwise  station  is  used  for  calculating  the  angle  of  attack,  and  it  is  given 
by 

V3/4(y)  =  V  +  S(u,s)x3/ 4  (7) 

where  x3/4  =  \{xa  —  0.5)c  0  0]T.  Let  [u3/4(y),  v3/4 (y),  w3/4{y)\  denote  the  components  of  V3/4  in 
the  local  station  frame,  and  V3/4  its  magnitude.  Then,  the  local  angle  of  attack  and  sideslip  can  be 
calculated  using 

tan  a(y)  =  —^~,  sin  j3(y)  =  (8) 

u3/4  V3/4 

5.4  Equations  of  Motion  of  the  Aircraft 

Let  mR  and  mR  denote  the  masses  of  the  right  wing  and  left  wing,  respectively.  Let  fhR  and  mR 
denote  the  masses  per  unit  length  of  the  right  wing  and  left  wing,  respectively.  Let  m  denote 
the  total  mass  of  the  aircraft,  and  let  tcg  denote  the  position  of  its  center  of  gravity.  Let  rs  = 
[xec,  0, 0]T  denote  location  of  the  center  of  gravity  with  respect  to  the  wing  twist,  and  let  y  = 
[xeccos0s(y),  y,  £  —  xecsin#s(y)]T  denote  the  position  of  the  center  of  gravity  of  a  wing  station  in 
the  wing  root  frame.  Let  u>s  =  [0  9S  0]T  denote  the  angular  velocity  of  a  given  wing  station  due  to 
twisting. 
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The  translational  equations  of  motion  are  given  by 

m  (ub  +  S(ljb)ub  +  S(cjB)rcG  +  S2(uB)rcG  +  S(ws)rcG)  + 

fb/ 2 

inn  /  ((S(u)r  +  u>b)S(u>r)  +  S(u>R))TBRy  +  S(ujR)TBRUf)  dy 


Jo 


+mL  f  {(S{uL  +  uB)S(uL)  +  S(uL))  TBrY  +  S(uL)TBRuf)  dy 
J -6/2 


fW  2 


+rriR  /  (TBRUf(y)  +  S(wb  +  Ur  +  TBRU>s)TBR(uf  +  S(us)rs)  +  TBrS(us)ys)  dy 


Jo 


+rhL  [  (TBRUf(y)  +  S(wB  +  wl  +  TBRWs)TBR{uf  +  S(wa)rs)  +  TBR{S(us)rs)  dy 
J  -6/2 


(9) 


[XB  Yb  ZBf 


where  [XB,  Yb  ,%b\  is  the  net  force  acting  on  the  aircraft  (aerodynamic  plus  gravitational),  with 
components  in  the  body  frame.  An  expression  for  the  net  force  is  given  Sec.  5.7.  The  position 
vector  of  the  center  of  gravity  (CG)  is  given  by 


1  (  fb/2  [°  \ 

r cg  =  ~  rnR  /  (uf  +  S(ujR)TBRy)dy  +  mL  (u/  +  S(u>L)TBRy)dy 
m  \  Jo  J-b/2  ) 


(10) 


For  highly  flexible  or  rapidly  flapping  wings,  the  dynamics  of  the  CG  serve  to  couple  the  transla¬ 
tional  and  rotational  dynamics  tightly.  The  CG  location  can  be  changed  using  an  actuated  mass, 
such  as  the  bob  weight  in  Doman,  Oppenheimer,  and  Sigthorsson  [23],  for  controlling  the  aircraft 
attitude. 

The  moment  of  inertia  of  the  right  wing  is  given  by 

fb/  2 

J r  =  —  S  (TRRy  +  TRRX.)dm  (11) 

Jo 


and  Jl  is  defined  similarly.  Let  pwjt  and  pWtL  denote  the  densities  of  the  right  and  left  wing, 
respectively,  and  J  is  the  total  moment  of  inertia  of  the  aircraft.  It  has  been  assumed  that  the 
moment  of  inertia  of  the  wing  is  constant  in  magnitude,  i.e.,  the  effect  of  wing  bending  and  twist 
on  the  net  moment  of  inertia  of  the  aircraft  is  ignored.  Subject  to  this  assumption,  the  following 
dynamic  equation  for  rotational  motion  is  obtained: 


J  Cjb  +  S(ujb)3ojb  +  m5(ufl)S(rcG)uB  +  mS(rcG)uB  +  S(rcG)uB  +  J  rojr  +  J  lojl 
+S(ujb)(JrUr  +  J  l^l)  +  (S(ljr)J  r  —  J  rS(ujr))(ujb  +  +  (S(ur)Jr  —  JrS(ur))(ujb  +  Ul) 

rb/2 


- rriR 


-mL 


J  ^ S(uB  +  uR)S(TBRy)TBR{uf  +  S(  Ts)us)  +  5'(Tsfiy)TSfl(u/  +  5(r3)ws)^  dy 
J  ^S(uB  +  uL)S(TBRy)TBR(uf  +  S( rs)us)  +  S(TBRy)TBR(uf  +  S{ra)djs)^j  dy 


rb/2 


+  /  {Pw,r{TbrJsus  +  S(ub  +  ur)TBrJsUs)  —  rhRS(TBRUf)TBRS(rs)u>s)dy 

Jo 

f° 

+  /  {Pw,l{T/br3sUs  +  S(ub  +  ul)TbrJsus))  —  fhLS(TBRUf)TBRS(rs)<jjs)dy  =  [L  M  TV]1 
J  ~ 6/2 


(12) 


where  J s(y)  =  —  f  fs  S2(x)dA  denotes  the  second  moment  of  area  matrix  of  a  cross  section  of 
the  wing.  An  expression  for  the  net  moment  ([ L  M  IV])  is  given  in  Sec.  5.7.  Note  that  if  the 
terms  arising  from  flexibility  are  ignored  along  with  the  wing  root  angular  velocity,  then,  with  the 
additional  assumption  that  yqg  =  0,  Euler’s  equations  are  recovered  as  one  would  expect.  The 
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equations  of  motion  derived  in  this  equation  incorporate  wing  rotation  (see  Eq.  (5)  which  expresses 
u>R  in  terms  of  the  flapping  rates)  and  therefore,  this  model  can  be  used  for  a  study  of  flexible 
flapping  wings  as  well.  These  equations  of  motion  can  be  specialized  for  a  rigid  aircraft,  as  described 
in  [64], 

5.5  Structural  Dynamics 

The  bending  and  twisting  elastic  equations  of  motion  for  the  right  wing  are  given  by 

inn  -rhnxec  1  [  [v]3  1  [  v{EIb£")"  +  (EIb£")"  -  T£"  1  =  r  Fs  S  1 

- ihRXec  Ip  \  fi  +  -v  (GJejy  -  (Gj0'a)'  [  Ms,2 

where 

D  =  ll)s  +  T  sb(ub  +  ur), 

and 

V  =  TSBuB  +  \if  +  S(ujb  +  Ur)(TsBub  +  U/)  +  T  sbS{u>b  +  uR)TBRUf 
+T5S  (S(uB  +  uR)  +  S2(ujb  +  uR))  TBRy 
£  =  £/  -  y$R 

Uf  =  [0  (U/]T,  (15) 

Remark:  The  displacement  £  should  be  viewed  as  comprised  of  the  deformation  £j,  and  a  rigid 
component,  ySR,  i.e.,  £  =  £/(y)  — ydR ,  with  £j(0)  =  0.  This  perspective  is  helpful  from  the  point  of 
view  of  practical  implementation  of  boundary  control  schemes.  Likewise,  one  may  consider  6S  as  the 
sum  of  flexible  and  rigid  twist  contributions  (denoted  by  0R).  instead  of  a  pure  deformation.  Then, 
the  wing  may  be  viewed  as  being  clamped  at  the  root,  with  9S(0)  =  6 R  +  0  (zero  deformation  at  the 
root).  This  decomposition  of  £  and  9S  changes  neither  the  governing  equations  nor  the  boundary 
conditions,  because  the  rigid  terms  do  not  affect  the  stiffness  and  damping  terms,  while  they  are 
already  incorporated  into  the  accelerations  and  the  right  hand  side. 

Note  that  Ip  =  pw Js(2,2)  and  =  Js(l,l),  where  pw  denotes  the  density  of  the  wing.  Fur¬ 
thermore,  F  S;3  =  FSi3(a,  d,  Fqo,  Uf,  9,  0)  is  the  total  force  acting  in  the  local  z  direction  (hence 

the  subscripts  ‘s’  and  ‘3’),  while  MS]2  —  MSi2(<T  oc,  Voo,  u/,  9, 0)  is  the  local  pitching  moment. 

The  arguments  of  F  and  M  listed  here  are  by  no  means  exhaustive;  rather,  they  are  the  primary 

contributors.  The  term  [V]3  denotes  the  ^-component  of  the  local  acceleration,  and  [u)\r2  is  the 
y-component  of  the  local  angular  acceleration.  Expressions  for  the  net  force  and  moment  are  given 
in  Sec.  5.7.  The  Kelvin- Voigt  damping  coefficient  is  obtained  by  scaling  Elb  and  GJ  by  a  factor  of 
r]  in  the  bending  and  twist  equations,  respectively. 

Remark:  The  scaling  term  7/  will  not  be  equal  for  both  cases,  viz.,  bending  and  twist,  in  the 
most  general  case.  Furthermore,  it  is  common  among  structural  dynamicists  to  model  the  damping 
coefficient  as  a  linear  combination  of  the  mass  (or  moment  of  inertia)  and  stiffness. 

Remark:  The  linear  model  presented  here  can  be  readily  replaced  by  a  nonlinear  model  in  the 
proposed  coordinates  to  match  the  requirements  of  the  problem  at  hand. 

The  boundary  conditions  are  given  by  the  following  expressions. 

•  At  the  wing  root:  £  =  0,  while  £;  and  9S  can  be  set  arbitrarily  (within  admissible  limits)  as 
the  dihedral  angle  and  the  twist,  respectively,  at  the  wing  root. 


(13) 

(14) 
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•  At  the  wing  tip:  =  0,  (. El {,£")'  —  T£'  =  0  and  6's  =  0  (i.e.,  free  end  boundary  conditions). 

Remark:  If  the  tension  T  is  spatially  varying,  i.e.,  T  =  T(y),  then  an  additional  term,  T'Q,  needs 
to  be  added  alongside  T £"  in  Eq.  (13). 

Boundary  conditions  at  the  wing  root,  in  particular  6S{ 0)  and  £'(0),  can  be  controlled  actively 
via  dedicated  actuators  for  stabilizing  an  unstable  wing  or  for  ensuring  that  the  net  force  on  the 
wing  or  its  components  achieve  the  desired  value  for  specific  maneuvers  [63]. 

5.6  Fuselage  Kinematics 

The  fuselage  attitude  is  described  by  the  Euler  angles  p,  6  and  p.  The  kinematic  equations  are 
given  by 

p  =  p  +  q  tan  9  sin  (j>  +  r  tan  9  cos  p, 

9  =  qcos(j)  —  r  sinpt  (16) 

p  =  (q  sin  <p  +  r  cos  p)  sec  9 

The  equations  which  relate  the  position  of  the  aircraft  to  its  translational  velocity  are  essentially 
decoupled  from  the  flight  dynamics,  and  are  given  by 

X  =  Vgn  cos  7  cos  y 

Y  =  tgn  cos  7  sin  y  (17) 

Z  =  —  Vgn  sin  7 

where  Vgn  is  the  ground  speed  of  the  aircraft  (obtained  by  subtracting  the  velocity  of  the  wind  from 
that  of  the  aircraft).  The  flight  path  angle  (7)  and  the  wind  axis  heading  angle  (y)  in  Eq.  (17)  are 
defined  as  follows: 


sin  7  =  cos  a  cos  /3  sin  9  —  sin  (3  sin  p  cos  9  —  sin  a  cos  (3  cos  p  cos  9 
sin  y  cos  7  =  cosacos(3cos9sinp  +  sin(3(smpsm9smp  +  cospcosp) 
+  sin  a  cos  (3{ cos  p  sin  9  sin  p  —  sin  p  cos  p) 

The  turn  rate  is  given  by  oj  =  y.  If  9  =  p  =  $  =  0,  it  follows  that 

oj  =  p  =  (q  sin  p  +  r  cos  p)  sec  9 


(18) 


(19) 


5.7  Forces 

The  net  force  on  the  aircraft  consists  primarily  of  contributions  from  aerodynamic  and  gravitational 
forces.  For  numerical  analysis,  the  aerodynamic  forces  and  moments  are  calculated  using  strip 
theory.  As  a  method,  strip  theory  is  used  for  aircraft  aeroelastic  simulations  [108]  and  routinely  for 
blade  element  theory  in  the  rotorcraft  field  [37].  Strip  theory  approaches  have  also  been  applied  to 
wings  in  a  trailing  vortex  flow  and  aircraft  spin  prediction  (see  references  [50,  69,  59,  58]  and  others 
cited  therein). It  seems  that  only  recently  has  the  general  strip  theory  approach  been  applied  in 
realtime  simulation  for  fixed-wing  force  and  moment  calculations  [94,  39,  38,  83,  84],  Strip  theory 
methods  have  also  been  applied  to  flapping  wing  aircraft  [20,  46,  14]. 
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The  wing  is  divided  into  chord-wise  strips.  The  lift,  drag,  and  the  quarter  chord  aerodynamic 
moment  at  each  strip  can  be  computed  by  using  a  suitable  aerodynamic  model,  and  these  can  be 
summed  over  the  entire  wing  to  yield  the  net  aerodynamic  force  and  moment.  A  similar  calcula¬ 
tion  is  performed  for  the  horizontal  tail  and  added  to  the  wing  contributions.  The  aerodynamic 
contributions  of  the  fuselage  are  ignored  with  the  understanding  that  they  can  be  added  readily  to 
the  model. 

Since  the  model  developed  is  intended  to  be  as  generic  as  possible,  the  model  proposed  by 
Goman  and  Khrabrov  [29]  is  presented  in  this  section  as  a  candidate  model  for  computing  the  lift 
and  the  quarter  chord  moment  while  drag  is  estimated  assuming  the  classic  drag  polar.  In  the 
authors’  estimate,  Goman  and  Khrabrov’s  model  offers  at  least  two  advantages  over  the  existing 
models  (such  as  Theodorsen  [97]  or  Peters  [68]).  First,  the  model  is  cast  in  the  form  of  a  single 
ordinary  differential  equation  (ODE)  and  two  algebraic  equations,  one  each  for  lift  and  the  quarter 
chord  pitching  moment.  The  state  variable  for  the  ODE  corresponds,  physically,  to  the  chordwise 
location  of  flow  separation  on  the  airfoil.  Therefore,  the  model  is  quite  easy  to  implement  as  part  of 
a  numerical  routine.  Second,  the  model  is  inherently  nonlinear  and  applicable  to  post-stall  flight. 

The  following  equation  describes  the  movement  of  the  separation  point  for  unsteady  flow  con¬ 
ditions 

T\v  +  v  =  v0(a  —  T2d)  (20) 

where  v  denotes  the  positon  of  the  separation  point,  n  is  the  relaxation  time  constant,  and  r2 
captures  the  time  delay  effects  due  to  the  flow,  while  vq  is  an  expression  for  the  nominal  position 
of  the  separation  point.  These  three  parameters  need  to  be  identified  experimentally  or  using  CFD 
for  the  particular  airfoil  under  consideration.  The  coefficients  of  lift  and  quarter-chord  moment  are 
then  given  by 


Ci  =  ^  sin(a(l  +  u  +  2y/u)) 

5  +  5^-6071 

16  J  v  ; 

The  lift  force  and  the  quarter  chord  moment  per  unit  span  are  then  given  by 
L{y)  =  0.5 pV(y)2cCl  +  ^pc2  (£ +  V^a  -  (xa  -  0.25)caj 

M(y)  =  0.5 PV{y)2c2C*mac  +  ^pc2  (v^  +  {Xa  ~  -  c2  +  (*„  -  0.25)2)  (22) 

where  a  is  the  local  angle  of  attack,  and  p  denotes  the  density  of  air.  Furthermore,  V  =  ||V||  is 
the  local  wind  speed  with  V  defined  in  Eq.  (6),  and  I^o  is  the  freestream  speed  of  the  aircraft  given 
by  Voo  =  || ||-  The  last  term  of  each  expression  was  added  to  Goman’s  original  model  [29]  and 
corresponds  to  the  apparent  mass  effect  [20]. 

There  is,  unfortunately,  no  simple  expression  for  the  sectional  drag  coefficient.  The  sectional 
drag  coefficient  can  be  written  as 


C*mac  =  \  sm(a(l  +  v  +  2^)) 


Cd 


0.664  1  2 

VI =e+^A~R  l' 


(23) 


where  Ci 


L(y ) 

0.5  pV{y)2c 


Ar  is  the  aspect  ratio  of  the  wing,  Re  denotes  the  chordwise  Reynolds 
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number,  and  e  is  Oswald’s  efficiency  factor.  The  skin  friction  term  [45]  assumes  laminar  flow 
over  the  wing  and  may  need  to  be  replaced  with  a  different  approximation  (see  DeLaurier  [20]  for 
instance).  The  drag  model  is  quasi-steady  in  nature  so  that  dynamic  stall  effects  are  not  included. 
A  refined  model  for  calculating  drag,  incorporating  dynamic  stall,  may  be  found  in  DeLaurier  [20]. 
The  local  aerodynamic  force  on  each  wing  can  be  written  in  the  body  axis  system 


'  XA(y)  ' 

L(y)  sin  a(y)  -  D(y)  cos  a(y) 

YA(y) 

=  T  BS 

0 

zA(y) 

_  ~L(y)  cos  a(y)  —  D(y)  sin  a(y)  _ 

(24) 


Note  that  T bs  is  used  instead  of  T br  and  this  is  the  most  important  source  of  the  difference 
between  the  resultant  of  the  forces  and  moments  on  a  flexible  wing  vis-a-vis  a  rigid  wing.  The 
components  of  the  gravitational  force  are  given  by 


Xg  =  —mg  sin  9,  Yg  =  mg  cos  9  sin  </>,  Zg  =  mg  cos  9  cos  4> 

and  the  corresponding  moment  is  given  by  S(rcg)[Xg  Yg  Zg]T . 

The  net  aerodynamic  force  on  two  wings  is  given  by 


'  xB  ' 

Yb 

nb/2 

'  XA(y)  ' 
YA(y) 

dy  + 

j-b/2 

'  XA(y)  ' 
YA(y) 

Zb 

Jo 

wing 

_  ZA{y)  _ 

right 

Jo 

ZA(y ) 

The  net  aerodynamic  moment  due  to  the  two  wings  is  given  by 


L 

rb/2 

'  XA{y)  ' 

rb/2 

'  XA(y)  - 

M 

= 

/  -S’(y) 

YA(y) 

dy  + 

/  -S’(y) 

YA(y) 

N 

wing 

Jo 

_  ZA(y)  _ 

right 

Jo 

ZA(y ) 

(25) 


(26) 


(27) 


A  similar  calculation  can  be  performed  for  the  horizontal  tail.  The  net  force  and  moment  on  the 
aircraft  themselves  are  the  sum  of  the  contributions  from  the  wing,  the  horizontal  tail  and  gravity: 


XB  ' 

'  vB  ' 

'  XB  ' 

r  ^  1 

Yb 

= 

Yb 

+ 

Yb 

+ 

Yg 

Zb 

Zb 

wing 

Zb 

tail 

L  Zg  \ 

L 

"  L 

L 

M 

= 

M 

+ 

M 

+  S(  rcg) 

N 

N 

wing 

N 

tail 

L  Zg  J 

This  completes  the  formulation  of  the  equations  of  motion. 


5.8  Aerodynamic  Model  for  Trim  Calculations 

The  lift  and  drag  coefficients  of  the  wing  and  tail  airfoils,  adapted  from  values  determined  experi¬ 
mentally  [100]  for  the  Vapor  airplane  itself  at  a  Reynolds  number  of  20000,  are  given  by 

CL  =  0.28295  +  2.00417a,  CD  =  0.0346  +  0.3438C£,  (29) 
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Angle  of  Attack  (deg) 


Figure  6:  Experimentally  obtained  aerodynamic  data  [100]. 


where  a  is  measured  in  radians.  Using  thin  airfoil  theory  [49],  it  was  determined  that  Cmac  = 
—0.1311.  The  actual  experimental  plot  has  been  shown  in  figure  6.  The  Cl  and  Cp  expressions  in 
equation  (29)  are  obtained  by  averaging  only  over  the  red  points  in  figure  6.  The  rest  of  the  points 
in  figure  6  (marked  in  black)  represent  data  collected  at  high  values  of  pitch  rate  and  a,  and  are  not 
relevant  to  the  discussion.  During  the  experiments,  for  a  >  25  deg,  a  was  seen  to  be  substantial, 
and  therefore,  the  coefficients  in  equation  (29)  are  reliable  only  up  to  a  =  25  deg. 

The  aircraft  weighs  12  grams,  including  a  ballast  mass  added  to  the  nose  of  the  aircraft  for 
placing  the  CG  around  half-wing-chord  under  nominal  conditions;  i.e.,  when  the  wing  dihedral 
and  incidence  are  both  zero.  The  aircraft  is  29.7  cm  long  from  nose  to  tail,  and  under  nominal 
conditions,  the  distance  between  the  AC  and  the  CG  is  xac  =  3.6  cm.  The  horizontal  tail  is  located 

26.1  cm  behind  the  wing  root  AC.  The  limiting  value  of  the  horizontal  tail  deflection  is  assumed  to 
be  30  deg  in  both  directions.  The  limiting  value  of  the  wing  dihedral  is  assumed  to  be  60  deg  on 
either  side,  while  that  of  the  wing  incidence  is  15  deg. 

6  Performance  and  Stability  of  a  Rigid  Aircraft 

6.1  Comparison  With  the  Vertical  Tail 

Figure  7  illustrates  the  physics  underlying  the  use  of  wing  dihedral  as  a  control.  Increasing  the 
wing  dihedral  reduces  the  force  acting  in  the  body  z-direction,  and  generates  a  side  force.  The 
reduced  z-force  affects  the  aircraft  flight  path  angle  and  angle  of  attack,  and  hence  the  flight  speed. 
On  the  other  hand,  the  side  force  can  be  used  for  providing  the  centripetal  force  for  turning,  and 
as  a  source  of  the  yawing  moment.  In  particular,  if  the  CG  is  located  behind  the  line  of  action  of 
the  side  force,  then  a  positive  side  force  produces  a  positive  yawing  moment  and  vice-versa  (see 
figure  5  for  the  sign  conventions).  It  follows  that  a  positive  rolling  moment  (wherein  the  lift  on 
the  left  wing  is  higher  than  the  right  wing)  is  accompanied  by  a  positive  yawing  moment  if  the 
wings  have  a  positive  dihedral  deflection.  Consequently,  the  adverse  yaw  produced  due  to  rolling 
is  reduced. 

Assuming  a  linear  relation  between  the  lift  and  the  angle  of  attack,  the  yawing  moment  generated 
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Reduce  net  lift  on  aircraft 
without  altering  drag 


Excess  lift  on  one  wing  produces 
side  force  and  yawing  moment 


Larger  dihedral  on  one  side 
produces  side  force  and  yawing 
moment 


Figure  7:  Illustration  of  the  physics  underlying  the  use  of  dihedral  as  a  control.  The  dark  conspic¬ 
uous  dot  in  the  figures  is  the  aircraft  CG. 


by  a  dihedral  deflection  5  of  the  left  wing,  while  that  of  the  right  wing  is  zero,  is  given  by 

Nw  =  qooSw(lwCLaa  +  cCniac)5  (30) 

Clearly,  the  dihedral  is  more  effective  for  yaw  control  at  high  angles  of  attack.  Equation  (30)  also 
suggests  that  the  dihedral  is  better  than  the  vertical  tail  when  It  is  small.  The  ability  to  change 
wing  dihedral  is  built  into  birds  in  the  form  of  their  ability  to  flap  their  wings  for  propulsion. 
Hence,  no  additional  mechanisms  are  needed  for  yaw  control.  Ornithopters,  too,  can  benefit  from 
differential  dihedral-based  yaw  control  in  a  similar  manner. 

For  positively  cambered  wings,  Cmac  <  0.  Hence  the  second  term  on  the  right  hand  side  of  (30) 
is  negative  and  not  only  reduces  £,  but  could  also  render  it  negative.  This  leads  to  the  problem  of 
reversal  of  control  effectiveness  (measured  by  the  sign  of  Nw  in  equation  (30)). 

Figure  8  plots  the  sign  of  the  control  effectiveness,  i.e.,  sign  (  a(5^N5r))  ,  on  a  p  —  r  grid  for 

a  =  8.595  deg  (0.15  rad).  The  plot  clearly  shows  that  the  sign  of  the  control  effectiveness  depends 
strongly  on  the  angular  rates.  Moreover,  as  shown  in  [64],  it  is  uniformly  negative  for  a  <  6  deg 
and  uniformly  positive  for  a  >  12  deg.  The  sign  of  the  control  effectiveness  is  usually  assumed 
to  be  known  a  priori  while  designing  control  laws.  The  challenge  involved  in  designing  a  sound 
turning  flight  controller  is  captured  in  figure  8. 

6.2  Analytical  Approximations  to  Lateral-Directional  Stability 

The  lift  and  drag  forces  produced  by  the  wing  as  well  as  their  moments  about  the  origin  of  the  body 
frame  can  be  resolved  along  the  body  axes.  In  particular,  summing  the  body  axis  components  of 
the  net  moment  due  to  lift  and  drag  yield  the  net  rolling,  pitching  and  yawing  moments.  Stability 
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Figure  8:  Plot  showing  the  sign  of  the  control  effectiveness,  sign  ^ /\(8L-8R) )  as  a  function  of  roll 
rate  and  yaw  rate  at  a  ss  8.5  deg. 


of  aircraft  depends  primarily  on  the  three  aerodynamic  moments  and  their  derivatives  with  respect 
to  the  aircraft  angular  velocity,  angle  of  attack  and  sideslip. 


Table  1:  Stability  derivatives  for  a  tailless  aircraft  with  an  articulated  wing 


Derivative 

Symmetric  flight 
(6l  =  Sr  =  5) 

Turning  flight 

(SL  «  SR) 

Stability  condition 

Lf> 

Stable 

Stable  when  Sr  +  Sr  >  0 

L/3  <  0 

Lp 

Stable 

Stable 

Lp  <0 

Np 

Unstable 

Unstable 

Np  >  0 

Nr 

Stable 

(drag  reduces  stability) 

Unstable,  but  stable  when 
sign(p)  /  sign(hL  -  Sr) 

Nr  <  0 

Table  1  lists  some  stability  derivatives,  and  their  qualitative  properties  for  a  tailless  articulated 
wing  aircraft.  Based  on  the  results  in  table  1,  it  is  clear  that  the  aircraft  would  be  expected  to 
be  unstable  in  most  flight  regimes.  At  least  two  stability  derivatives  suggest  the  possibility  of 
stability  in  some  select  turn  regimes.  However,  in  rapid  turn  regimes,  the  flight  dynamics  are  far 
too  strongly  coupled  to  draw  reliable  conclusions  from  linear,  decoupled  analysis.  Furthermore, 
because  the  stability  derivatives  depend  strongly  on  the  wing  dihedral  angle,  which  is  in  turn  a 
function  of  the  aircraft  maneuver,  it  follows  that  stability  is  tied  very  closely  to  the  nature  of 
the  maneuver  being  executed.  The  observation  is  peculiar  to  aircraft  with  articulated  wings.  In 
conventional  aircraft  with  fixed  wings,  although  stability  derivatives  depend  on  aircraft  states,  they 
are  essentially  independent  of  the  control  surface  deflection. 

The  idea  of  using  wing  dihedral  for  control  is  particularly  useful  when  the  wings  are  flexible, 
because  flexible  wings  bend  and  twist  spontaneously  under  aerodynamic  loading.  The  dihedral 
angle  at  a  given  point  on  the  wing  is  equal  to  the  sum  of  the  slope  of  the  bending  displacement 
and  the  wing  slope  at  the  root.  Since  bending  and  twisting  are  coupled,  wing  twist  can  be  used  to 
bring  about  a  passive  proverse  change  in  the  wing  dihedral. 
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7  Analysis  of  the  Wing  and  Effective  Dihedral 

7.1  Effective  Dihedral 

The  Young’s  modulus  of  the  wing,  E,  may  be  considered  as  a  design  parameter.  In  order  to  exploit 
the  idea  of  using  the  wing  dihedral  for  yaw  control,  the  wing  dihedral  effect  itself  may  be  looked 
upon  as  a  design  driver  for  E. 

The  role  of  differential  (or  asymmetric)  dihedral  for  yaw  control  has  been  discussed  in  detail 
in  Ref.  [64],  The  dihedral  primarily  produces  a  side  force,  which  is  actually  a  component  of  the 
total  force  produced  by  the  wing  normal  to  its  local  plane.  Let  Ya  and  Za  denote  the  local  forces 
produced  by  the  wing  along  the  body  y  and  z  axes,  respectively.  Therefore,  one  may  define  a  term 
called  the  effective  dihedral,  <5eff ,  as  follows: 


Seff  =  tan 


-l 


(  Jo/2YA(y)dy\ 
\lS/2  ZA(y)dyJ 


(31) 


This  notion  of  effective  dihedral  is  different  from,  and  arguably  more  general  than,  that  of  Rodden 
[77]  who  derived  expressions  for  the  increments,  arising  from  the  wing  bending,  in  the  rolling 
moment  derivatives.  The  notion  of  effective  dihedral  is  particularly  useful  for  wing  design  from 
the  point  of  view  of  elasticity.  The  Young’s  modulus,  E,  could  be  chosen  to  ensure  that  the 
wing  produces  a  sufficient  effective  dihedral  effect  with  reasonable  actuator  forces.  The  effective 
dihedral  depends  on  the  boundary  conditions  to  which  the  wing  is  subjected  whereas  the  boundary 
conditions  themselves  depend  on  the  location  and  type  of  actuators.  For  a  rigid  wing,  the  effective 
dihedral  and  the  actual  dihedral  are  equal. 


(a)  Effective  dihedral  when  E  —  5  MPa. 


(b)  Effective  dihedral  when  E  —  50  MPa. 


Figure  9:  Effective  dihedral  as  a  function  of  the  dihedral  angle  at  the  wing  root  for  two  different 
values  of  the  Young’s  modulus.  Each  plot  shows  the  effective  dihedral  for  three  values  of  wing  tip 
twist  (6):  0,  0.1  rad  and  0.2  rad.  This  plot  was  obtained  for  V  =  2.5  m/s  and  a  =  10  deg. 

Figures  9(a)  and  9(b)  show  the  effective  dihedral  as  a  function  of  the  wing  dihedral  angle  at 
the  root.  The  effective  dihedral,  as  expected,  is  much  higher  for  E  =  5  MPa  as  compared  to 
E  =  50  MPa.  In  the  former  case,  the  wing  bending  is  large  enough  so  that  the  flexibility  provides  a 
substantial  increase  in  the  wing  dihedral  effect.  This  suggests  that  for  the  particular  wing  geometry 
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considered  here,  a  material  with  a  Young’s  modulus  of  E  ~  0(1)  MPa  should  be  chosen  in  order 
to  obtain  a  significant  dihedral  effect.  This  conclusion  depends  on  other  chosen  parameters  and 
hence,  such  analysis  should  be  performed  on  a  case-by-case  basis.  Furthermore,  it  is  important  to 
note  that  the  effective  dihedral  depends  on  the  trim  condition  under  consideration. 

The  effective  dihedral  is  useful  in  another  way.  It  forms  the  basis  to  extend  the  stability  analysis 
for  a  rigid  aircraft  to  the  case  of  flexible  wings.  In  Ref.  [64],  for  example,  analytical  expressions 
for  the  traditional  lateral  stability  derivatives  were  obtained  for  a  rigid  aircraft  and  the  stability 
of  lateral-directional  modes  was  examined  for  various  values  of  the  wing  dihedral.  Those  results 
would  be  applicable  to  a  flexible-winged  aircraft  when  the  effective  dihedral  angle  of  the  wing  is 
matched  to  the  dihedral  angle  of  a  rigid  wing.  This  is  valid  regardless  of  the  deformation  profile 
of  the  wing.  For  the  aircraft  model  considered  here,  it  suggests  that  the  motion  stability  would  be 
similar  to  that  of  the  rigid  aircraft  when  E  >  0(10)  MPa. 


7.2  Bending  and  Twist  Natural  Frequencies 

Traditionally,  natural  frequencies  of  lifting  surfaces  are  defined  in  terms  of  inertia  and  elastic  stiff¬ 
ness.  However,  unsteady  aerodynamic  lift  and  moment  relations  contain  terms  which  mathemati¬ 
cally  play  the  same  role  as  stiffness,  damping,  and  inertia  in  the  governing  relations.  Consequently, 
another  set  of  natural  pseudo  frequencies  can  be  defined  which  include  these  aerodynamic  contri¬ 
butions. 

Consider  the  case  where  0'(b/ 2)  =  £"(b/ 2)  =  £'"(b/2)  =  0.  If  ujq  and  denote  the  frequencies 
of  the  first  (decoupled)  twisting  and  bending  modes,  respectively,  then  it  can  be  shown  that  [6] 


2  7 r2  GJ  M  2  12.36  El ), 

Ue  =  AL2!^  ~  V  ^  =  ~1A 


(32) 


where  M  denotes 


9Ms,2 

do 


(linearized  twisting  moment).  In  order  to  estimate  the  extent  of  time-scale 


separation,  the  ratio  is  of  interest.  Time-scale  separation  is  a  property  wherein  the  dynamics 

consist  of  two  sets  of  modes,  one  of  which  is  significantly  (at  least  an  order  of  magnitude)  faster 
than  the  other  mode.  The  stability  of  each  mode  can  then  be  analyzed  independently,  with  other 
mode  contributing  a  constant  term  whose  value  is  a  function  of  the  mode  being  analyed  [1]. 

It  can  be  shown  that  the  time  scale  separation  between  bending  and  twisting  depends  primarily 
on  the  aspect  ratio,  the  flight  speed  and  the  geometry  of  the  wing  cross  section.  The  ratio  w|/w| 
is  given  by  [62,  63] 


_  G  L2  _  3  L2 

ooj  ^ e  c2  2(1  +  Up)  c2 


(33) 


where  vp  is  Poisson’s  ratio.  The  ratio  1.5/ (1  +  Vp)  ~  1.  Therefore, 

—  ps  0(L/c )  (the  aspect  ratio  of  the  wing) 

Therefore,  the  twist  dynamics  are  faster  than  the  bending  dynamics.  The  time  scale  separation 
reduces  with  decreasing  E  and  increasing  V,  as  the  influence  of  the  aerodynamic  terms  increasingly 
dominates  the  contribution  from  elasticity.  The  time  scale  separation  increases  with  increasing 
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aspect  ratio.  Although  this  time  scale  separation  cannot  be  used  to  draw  any  inference  about 
the  susceptibility  of  the  wing  to  flutter,  it  can  be  used  as  the  basis  for  designing  independent 
controllers  for  controlling  wing  bending  and  torsion.  Time  scale  separation  is  used  quite  commonly 
in  flight  control  design  in  a  similar  fashion  [89,  103,  73,  35].  Incidentally,  equation  (33)  suggests 
that  significant  time  scale  separation  is  to  be  expected  in  slow  aircraft  with  large  aspect  ratio  wings, 
such  as  gliders  and  high  altitude,  long  endurance  aircraft  like  the  NASA  HALE  [71,  72], 

8  Bifurcation  Analysis  of  Turning  Flight 

The  performance  and  stability  of  an  MAV  equipped  with  flexible  wings  (E  =  5  MPa)  in  steady 
turning  flight  is  analysed  in  a  manner  similar  to  that  described  for  a  rigid  aircraft  in  Ref.  [64],  A 
similar  analysis  could  be  repeated  for  other  maneuvers  of  interest.  Insofar  as  turning  is  concerned, 
wing  flexibility  may  have  one  or  more  of  several  possible  consequences. 

1.  The  overall  turn  rate  may  improve  because  of  the  additional  dihedral  generated  by  the  flexible 
wings. 

2.  Alternately,  for  a  given  turn  rate,  the  dihedral  angles  required  at  the  wing  root  would  be 
reduced. 

3.  When  the  sideslip  is  not  deliberately  regulated,  it  would  be  reduced  due  to  the  enhanced 
dihedral  effect. 

It  turns  out  that  flexibility  does  result  in  a  net  improvement  in  the  turn  rate  of  the  aircraft,  but 
only  when  wing  incidence  angle  at  the  root  (or  wing  twist  in  general)  is  used  actively.  There  is  a 
significant  reduction  in  the  sideslip  when  the  wings  are  locked  in  a  symmetric  dihedral  configuration. 
However,  when  the  dihedral  angles  alone  are  used  for  turns,  the  maximum  achievable  turn  rate 
does  not  improve  vis-a-vis  a  rigid  aircraft.  Furthermore,  the  magnitude  of  the  commanded  dihedral 
deflections  required  for  a  given  turn  rate  is  reduced  in  comparison  to  an  aircraft  with  rigid  wings. 
The  complete  set  of  trim  analysis  may  be  found  in  [64,  62],  while  some  sample  cases  are  presented 
here. 

8.1  Reduction  in  Sideslip  (Variable  0r\  Or  =  —6l  =  —0a;  Sr  =  Sr) 

A  turn  is  usually  initiated  by  rolling  the  aircraft  to  the  appropriate  bank  angle  and  followed  by  pro¬ 
viding  the  appropriate  yaw  rate  and  pitch  rate.  When  the  flexible  wings  are  twisted  asymmetrically, 
the  resultant  roll  rate  causes  a  build-up  in  yaw  rate  due  to  the  dihedral  effect.  However,  if  the  wings 
are  locked  in  a  symmetric  dihedral  configuration,  the  resultant  turn  is  accompanied  by  a  sideslip 
which  increases  with  increasing  turn  (roll)  rate.  This  phenomenon  has  been  captured  in  Fig.  10 
where  the  dihedral  angle  at  the  root  was  set  to  29  deg  (0.5  rad)  for  both  wings.  The  equilibrium 
points  are  marked  with  a  red  asterisk  A’,  indicating  that  they  are  unstable  with  a  single  positive 
real  eigenvalue.  For  a  rigid  wing,  the  sideslip  remains  less  than  5  deg  until  the  turn  rate  builds  up 
to  35  deg/s  (compared  with  nearly  70  deg/s  for  flexible  wings).  Thereafter,  the  aerodynamic  data 
used  here  is  insufficient  to  provide  accurate  trim  results.  In  general,  though,  the  sideslip  increases 
with  increasing  turn  rate  for  an  aircraft  with  a  rigid  wing.  On  the  other  hand,  when  the  wings  are 
flexible,  the  turn  rate  increases  sharply  with  increasing  wing  twist  and  furthermore,  the  sideslip 
peaks  at  just  over  10  deg  and  drops  thereafter  due  to  the  increasing  effective  dihedral  angle.  With 
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(a)  Turn  rate  for  the  rigid  wing  air¬ 
craft 


(c)  Turn  rate  for  the  flexible  wing  air¬ 
craft 
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(b)  Sideslip  for  the  rigid  wing  aircraft 


(d)  Sideslip  for  the  flexible  wing  air¬ 
craft 


Figure  10:  A  comparison  of  the  sideslip  and  turn  rate  as  functions  of  anti-symmetric  wing  twist 
for  otherwise  identical  airframes  equipped  with  rigid  and  flexible  wings.  The  wings  have  a  Young’s 
modulus  of  5  MPa.  The  equilibria  are  marked  with  a  red  asterisk  to  denote  that  the  Jacobian 
has  a  single  positive  real  eigenvalue.  In  both  cases,  the  dihedral  angle  at  both  wing  roots  was 
set  to  25  deg.  The  flight  speed  was  set  to  2.8  m/s  and  the  elevator  was  fixed  at  —11  deg  and 
$L  =  5r  =  29  deg  (0.5  rad). 


aerodynamic  data  that  is  accurate  for  larger  values  of  sideslip,  the  value  of  sideslip  at  the  peak 
is  liable  to  shift  from  that  obtained  with  the  present  model.  However,  the  peak  itself  occurs  due 
to  a  favourable  yawing  moment  which  comes  with  an  increasing  wing  dihedral.  Therefore,  a  peak 
would  be  expected  even  with  improved  aerodynamic  data,  unless  adverse  yawing  moment  from  the 
fuselage  causes  the  sideslip  to  keep  increasing  with  the  turn  rate. 

It  is  of  interest  to  note  that  the  topology  of  the  equilibrium  surface  depends  strongly  on  the 
wing  dihedral.  If  the  root  dihedral  angles  are  set  to  zero,  a  qualitatively  different  picture  emerges 
as  shown  in  [62].  Therefore,  the  open-loop  behavior  needs  to  be  understood  thoroughly  before  a 
turning  controller  is  designed. 
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8.2  Coordinated  Turn  (6jj  —  Or  —  0;  Sl,  Sr  variable) 

Figure  11  compares  the  turning  performance  an  aircraft  equipped  with  rigid  wings  with  that  of  one 
equipped  with  flexible  wings  with  Young’s  modulus  E  =  5  MPa  when  the  sideslip  is  required  to  be 
regulated  to  {3  =  0.  The  twist  angle  at  each  wing  root  is  set  to  zero,  i.e.,  Or  =  0l  =  0.  It  is  clear 
that  there  is  no  appreciable  increase  in  the  maximum  achievable  turn  rate.  However,  a  noticeably 
smaller  dihedral  deflection  is  required  at  the  wing  root  for  a  given  turn  rate  when  the  wings  are 
flexible,  as  expected.  The  stability  characteristics  seen  for  the  two  sets  of  aircraft  are  identical.  The 


(a)  Turn  rate  for  the  rigid  wing  aircraft 


(c)  Turn  rate  for  the  flexible  wing  aircraft 


(b)  Right  wing  dihedral  angle  for  the  rigid 
wing  aircraft 


(d)  Right  wing  dihedral  angle  for  the  flexible 
wing  aircraft 


Figure  11:  A  comparison  of  the  turn  rate  as  a  function  of  the  left  wing  dihedral  angle,  and  the  right 
wing  dihedral  angle  required  to  maintain  zero  sideslip,  for  otherwise  identical  airframes  equipped 
with  rigid  and  flexible  wings.  In  both  cases,  the  elevator  deflection  was  fixed  at  —11  deg,  and  Or  = 
0l  =  0.  The  flexible  wings  have  a  Young’s  modulus  of  5  MPa.  The  Jacobian  of  equilibria  marked 
by  pink  dots  have  three  eigenvalues  with  positive  real  parts:  one  real  and  a  complex  conjugate  pair. 
The  flight  speed  and  angle  of  attack  are  within  the  range  of  validity  of  the  aerodynamic  data. 


28 


points  marked  A,  B,  C  and  D  are  all  Hopf  bifurcations.  Evidently,  none  of  the  computed  equilibria 
possess  inherent  stability. 

Remark:  It  was  seen  in  Sec.  8.1  that  the  turn  rate  improved  for  a  flexible  wing  MAV,  accompa¬ 
nied  by  a  reduced  sideslip.  On  the  other  hand,  in  the  present  section,  there  is  a  deterioration  in  the 
coordinated  turn  performance,  measured  by  the  maximum  turn  rate,  when  the  wings  are  flexible. 
This  can  be  explained  as  follows.  At  the  angle  of  attack  considered  here,  the  wing  twists  upward 
(i.e.,  the  leading  edge  goes  up)  so  that  the  net  angle  of  attack  on  the  wing  is  higher  than  in  the  rigid 
case.  Therefore,  for  a  given  tail  setting,  the  aircraft  flies  at  a  lower  flight  speed  to  maintain  trim  in 
pitch.  The  reduced  speed  leads  to  a  reduction  in  the  net  lift,  which,  in  turn,  reduces  the  amount  of 
centripetal  force  available  to  sustain  rapid  turns.  Another  point  worth  noting  is  that  the  maximum 
achievable  turn  rate  depends  on  the  maximum  achievable  yawing  moment.  The  yawing  moment  for 
a  given  wing  incidence  setting  reaches  a  maximum  when  the  wing  dihedral  angle  is  45  deg,  or  when 
the  effective  dihedral  of  a  flexible  wing  equals  45  deg.  This  sets  another  fundamental  limitation  on 
the  maximum  achievable  turn  rate,  and  one  that  arises  out  of  the  sole  use  the  wing  dihedral  for 
turning. 

8.3  Discussion 

The  results  presented  above  yield  some  interesting  design  pointers.  The  wing  flexibility  can  be 
reduced  significantly  (up  to  0(10)  MPa  for  the  wing  geometry  and  flight  speeds  considered  here) 
without  achieving  a  substantial  improvement  in  the  coordinated  turn  rate  or  any  measurable  change 
in  the  effective  dihedral  angle,  although  a  considerable  saving  in  the  wing  mass  can  be  achieved  in 
the  process.  The  motion  stability  (notwithstanding  the  structural  stability  of  the  wing)  will  not 
be  markedly  different  from  that  of  a  rigid  configuration.  One  interpretation  which  follows  is  that 
flexibility  offers  only  a  limited  improvement  in  the  performance,  notwithstanding  savings  on  the 
wing  mass.  Alternately,  a  complete  aeroelastic  analyis  can  be  bypassed  as  long  as  the  flutter  and 
divergence  speeds  are  “considerably  larger”  than  the  prescribed  flight  speeds  (see  Sec.  7.2). 

These  conclusions  are,  by  no  means,  universally  valid  but,  when  used  judiciously,  can  achieve 
considerable  savings  in  the  computational  effort  invested  in  the  design.  In  a  recent  paper,  Baghdadi, 
Lowenberg,  and  Isikveren  [3]  observed  that  the  open  loop  stability  characteristics  did  not  change 
markedly  between  the  rigid  and  flexible  configurations  considered  in  their  paper.  This  is  in  keeping 
with  the  observations  in  this  AFOSR  project.  Nevertheless,  a  control  law  designed  using  a  rigid 
model  yielded  markedly  different  closed  loop  stability  characteristics  when  the  time  constants  of 
the  rigid  and  flexible  modes  were  close  to  each  other.  On  similar  lines,  Merrett  and  Hilton  [53] 
demonstrated  that  flutter  (motion  instabilities)  can  arise  in  high-speed  aircraft  due  to  transient 
maneuvers  such  as  accelerations  or  rapid,  instantaneous  turns. 

9  Experiments,  Control  Design  and  Perching 

The  purpose  of  the  experiments  was  (a)  to  design  control  laws  for  an  articulated  wing  aircraft, 
and  (b)  to  implement  them  with  a  view  of  understanding  their  effectiveness  as  well  as  limitations. 
This  section  also  discusses  the  perching  maneuver,  which  is  one  of  the  most  important  maneuvers 
a  flapping-wing  aircraft  would  be  expected  to  execute  in  the  gliding  phase.  Further  details  on  the 
experiments  may  be  found  in  [26]. 

The  experiments  described  here  were  performed  on  a  test  MAV,  shown  in  Figure  12,  which  is  a 
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(a)  Symmetric  dihedral  40  degrees 


(b)  Asymmetric  configuration 


Figure  12:  Representative  configurations  showing  the  asymmetric  dihedral  capability  of  the  wings. 
The  foam  table  on  which  the  aircraft  is  resting  is  not  part  of  the  airframe. 


modified  version  of  the  commerically  manufactured  ParkZone  Ember  2.  The  aircraft  weighs  44  g, 
and  has  a  wing  span  of  42  cm  with  an  aspect  ratio  of  4.5.  The  outboard  sections  of  both  wings 
were  free  to  rotate  from  a  maximum  45  deg  dihedral  to  minimum  —15  deg  anhedral  for  a  total  arc 
range  of  60  deg.  The  actuators  for  wing  dihedral,  it  may  be  recalled,  are  controlled  independently 
on  both  wings  for  yaw  stability  and  control.  The  VICON  motion-capture  system,  consisting  of 
16  infrared  cameras,  was  used  to  collect  flight  data,  in  particular  the  aircraft  position  and  spatial 
orientation,  at  100  Hz. 

9.1  Control  Law  Design 

Control  law  design  for  the  MAV  is  described  in  this  section.  The  control  law  has  a  two-tier 
hierarchical  structure  based  on  time-scale  separation  [103]  which  occurs  naturally  between  the  fast 
rotational  dynamics  and  the  slow  translational  dynamics: 

•  The  innermost  loop  commands  the  elevator  and  the  asymmetric  components  of  the  wing 
dihedral. 

•  The  outer  loop  commands  the  angle  of  attack  and  turn  rate  to  be  tracked  by  the  inner  loop 
based  on  flight  speed  and  turn  rate.  The  turn  rate  and  the  flight  path  angle  are  computed 
based  on  position  measurements. 

A  schematic  of  the  controller  is  shown  in  Fig.  13. 

The  control  design  used  for  experiments  can  be  justified  using  the  dynamic  inversion  (DI) 
approach  presented  by  Hovakimyan,  Lavretsky  and  Sasane  [35].  In  [26],  we  argued  that  the  DI- 
based  controller  can  be  simplified  to  standard  PI(D)  controllers,  and  exact  gain  tuning  laws  are 
derived  in  the  process.  The  point  to  be  noted  here  is  that  PI  and  PID  controllers  can  be  used  for 
nonlinear  systems.  Singular  perturbation  theory  can  be  used  to  show  that  the  tracking  error  bound 
is  of  the  order  1/e  [35].  To  address  robustness  concerns,  a  filter  can  be  added  in  line  with  the  small 
gain  theorem.  This  is  done  in  adaptive  control  methodologies  where,  instead  of  substituting  for 
f(t,x,u )  with  x,  f(t,x,u)  is  predicted  online  using  adaptive  algorithms  [60]. 
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Figure  13:  Schematic  of  the  controller,  where  x  denotes  the  aircraft  heading. 


(a)  Wind  axis  angles,  a  and  /3  (b)  Heading  angle,  tp 

Figure  14:  Simulated  time  histories  of  the  aircraft  in  a  disturbance-free  flight.  A  12  deg  (0.2  rad 
jump  in  the  angle  of  attack,  a,  is  commanded.  The  resulting  disturbances  are  rejected  by  the 
control  law.) 

9.2  Simulations 

PID  controllers  designed  using  Dl-inspired  tuning  were  simulated.  The  time  histories  obtained 
in  two  different  environments  have  been  plotted  in  Figs  14  (no  external  disturbances)  and  15 
(persistent  periodic  disturbances).  In  both  cases,  the  controllers  performed  satisfactorily.  The 
angle  of  attack  was  kept  above  11  deg  to  ensure  the  yaw  control  effectiveness  of  the  dihedral  was 
uniformly  positive. 

The  purpose  of  the  simulations  was  to  demonstrate  a  general  control  design  technique.  However, 
in  the  course  of  experiments,  we  were  able  to  make  reasonable  estimates  of  the  open  loop  dynamics. 
This  allowed  us  to  tune  controllers  without  resorting  to  a  Dl-inspired  scheme. 

The  longitudinal  dynamics  of  the  aircraft  were  seen  to  be  stable,  but  poorly  damped  for  a  > 
8  deg.  Around  a  =  15  deg,  the  elevator  effectiveness  saturated  and  higher  angles  of  attack  were 
unattainable  under  steady  state  conditions.  The  open  loop  response  was  measured  to  have  a  time 
period  of  1  s.  The  observed  reduction  in  the  amplitude  of  oscillations  was  used  to  approximate  the 
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(a)  Wind  axis  angles,  a  and  ft 


(b)  Heading  angle,  ip 


Figure  15:  Simulated  time  histories  of  the  aircraft  in  a  persistent  periodic  lateral-directional  dis¬ 
turbance  field.  A  12  deg  (0.2  rad  jump  in  the  angle  of  attack,  a,  is  commanded.) 


damping  coefficient  to  0.046.  The  open  loop  dynamics  can  be  written  in  the  form 

a  +  0.62d  +  40a  =  —40 6e  +  5.6 


(34) 


Therefore,  an  essentially  derivative-integral  controller  was  designed  for  the  tailless  configuration: 

rt 


Se(t)  =  0.14  —  ac  +  kde  +  ki  /  eadt, 

Jo 


(35) 


where  the  offset  of  0.14  rad  was  added  based  on  the  measured  Se  —  a  trims.  The  gain  ki  was  similar 
to  that  for  the  configuration  with  a  vertical  tail,  while  kd  =  0.217  is  chosen  so  that  the  damping 
coefficient  is  approximately  equal  to  0.7. 

Based  on  experimental  observations,  it  was  estimated  that  the  open  loop  yaw-rate  dynamics 
are  of  the  form 

r  +  2 (tvf  +  o?r  =  C  ~  —0.1,  u  «  2vr  (36) 

for  a  <  8  deg.  Thereafter,  the  yaw  dynamics  are  unstable  and  oscillatory  in  nature.  In  order  to 
account  for  the  actuator  time  delay  of  0.2  s,  a  lead  compensator  L(s)  was  designed  given  by  L(s)  = 
8(s  +  4.5) 


4.5(s  +  8) ' 


The  role  of  dihedral  control  is  regulation,  and  it  suffices  use  a  derivative  controller 


12(s  +  4) 

for  damping  addition,  with  the  derivative  filter  given  by  D(s)  =  — -  —  — — .  The  commanded 
anti-symmetric  dihedral  deflection  is  given  by 


-'asym 


=  kdD(s)L(s)er(s) 


(37) 


9.3  Perching  Guidance  Loop 

The  outer  control  loop  is  designed  to  ensure  rapid  changes  in  the  flight  path  over  a  short  duration. 
The  flight  path  angle  is  controlled  in  discrete  time  so  that  a  symmetric  dihedral  angle  is  commanded 
every  0.2  s  (which  is  equal  to  the  dihedral  acutator  time  delay).  The  commanded  dihedral  angles 
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are  given  by 


8r  =  $L  =  ' i  /  2  + 


/(«) 


CL(a) 


/(a)  tan 7c’  Co  (a) 

where  yc  is  the  commanded  flight  path  angle  which  is,  in  turn,  given  by 

h  ,  z  —  zi 


7C  =  tan  l(h) 


1  +  0.28125/i2 


,  h  = 


V(x~  xi)2  +  {y  -  yi )2 


(38) 


(39) 


Here,  xi,  m  and  z;  are  the  coordinates  of  the  desired  landing  point  on  the  ground,  or  a  point  in  the 
air  where  a  perching  command  is  to  be  sent  to  the  aircraft.  This  simple  controller  was  seen  to  be 
quite  effective  over  the  short  duration  of  the  experiments. 


9.4  Angle  of  Attack  Control  for  Perching 


(a)  a,  9  and  7  (b)  Flight  speed 

Figure  16:  Experimental  results  showing  the  longitudinal  flight  parameters.  In  particular,  a  settles 
down  at  the  desired  value  within  2  s. 

Figure  16  shows  the  experimentally-measured  longitudinal  flight  parameters.  For  these  experi¬ 
ments,  the  wing  dihedral  was  not  controlled  actively  which  caused  the  aircraft  heading  to  deviate 
steadily  from  a  straight  line.  An  angle  of  attack  of  5  deg  was  commanded  while  the  flight  speed 
and  flight  path  angle  were  not  controlled.  The  controller  for  the  tailless  aircraft  yielded  similar 
characteristics  as  the  one  with  a  vertical  tail. 

9.5  Lateral-Directional  Control  for  Perching 

In  the  aircraft  with  a  vertical  tail,  local  lateral  stability  was  achieved  using  a  simple  PID  controller. 
However,  in  several  flight  tests,  the  roll  rate  was  seen  to  build  up  due  to  the  dihedral  effect  and, 
without  wing  twist  or  ailerons,  could  not  be  compensated.  This  led  to  a  divergent  lateral-directional 
behavior  despite  local  stability.  Figure  17  shows  the  time  histories  for  the  case  where  the  lateral 
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(a)  Sideslip  and  velocity  heading  (b)  Euler  angles 

Figure  17:  Experimental  results  showing  the  sideslip,  velocity  heading  and  the  Euler  angles  mea¬ 
sured  during  a  yaw  control  test  of  the  aircraft  with  a  vertical  tail.  Parameters  appear  to  be 
regulating  during  the  short  experiment 


dynamics  were  seen  to  be  stable.  A  zero  heading  angle  was  commanded.  The  heading  angle  as 
well  as  sideslip  converge  to  small  values.  However,  the  transient  response  does  not  vanish  within 
the  limited  flight  duration.  Nevertheless,  the  yaw  rate  slows  significantly  by  the  end  of  the  flight 
indicating  good  closed  loop  stability  characteristics.  Lateral  control  of  a  tailless  configuration  is 
under  experimental  investigation. 

9.6  Flight  Path  Control  for  Perching 

An  effective  flight  path  controller  is  necessary  for  a  successful  perching  maneuver.  The  aircraft  must 
be  able  to  track  the  desired  flight  path  in  order  to  arrive  at  a  spatial  target  with  an  acceptable 
flight  speed  and  height.  The  PID  controller  gains  on  the  angle  of  attack  controller  were  tuned  to 
ensure  sound  tracking  characteristics  across  a  range  of  flight  path  angles.  The  flight  path  angle 
itself,  as  explained  in  section  9.3,  is  controlled  using  the  wing  dihedral  angles. 

9.7  Perching 
9.7.1  Bioinspiration 

Figure  18  shows  some  snapshots  of  an  owl  executing  a  perching  maneuver,  extracted  from  a  reputed 
BBC  documentary  called  the  Life  of  Birds,  and  the  clip  was  processed  in  Matlab.  The  longitudinal 
flight  parameters,  the  flight  path  angle  7,  the  body  axis  pitch  angle  6,  and  the  angle  of  attack  a, 
were  extracted  by  making  two  assumptions:  (a)  the  depth  calibration  was  assumed  to  be  unchanged, 
and  (b)  the  local  ground  level  was  assumed  to  be  approximately  horizontal.  These  flight  parameters 
have  been  plotted  in  Fig.  19.  The  perching  maneuver  is  seen  to  consist  of  two  phases:  a  gliding 
phase  to  bring  the  aircraft  to  a  suitable  position  with  respect  to  the  landing  spot,  followed  by  a 
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Figure  18:  Snapshots  showing  an  owl  in  various  stages  of  a  perching  maneuver,  from  BBC’s  Life  of 
Birds.  The  video  was  processed  using  Matlab. 


Figure  19:  Angle  of  attack,  flight  path  angle  and  pitch  angle  measurements  from  the  BBC  video  of 
a  perching  owl.  The  maneuver  was  completed  at  t  =  1.5  s. 


rapid  pitch  up  (usually  to  a  post-stall  angle  of  attack)  which  leads  to  an  instantaneous  climb  and  a 
rapid  deceleration.  Interestingly  enough,  an  identical  profile  is  obtained  if  the  perching  maneuver 
is  designed  by  optimization  [17,  18]:  in  particular,  the  elevator  switched  between  two  values,  one 
corresponding  to  a  low-a  flight  and  the  other  being  the  saturation  value  of  the  elevator  deflection 
which  produced  the  desired  post-stall  angle  of  attack. 

9.7.2  Experimental  Demonstration  of  a  Perching  Maneuver 

In  conjunction  with  the  guidance  controller,  a  perching  maneuver  is  executed  as  follows.  An 
appropriate  altitude  is  chosen  such  that  a  perching  command  is  sent  when  the  aircraft  crosses  it. 
This  value  was  chosen  to  accommodate  the  actuation  time  delays  for  the  wing  dihedral  as  well 
as  the  elevator.  Until  this  point,  the  angle  of  attack  and  flight  path  angle  controllers  described 
in  Sec.  9.1  were  used  actively.  Once  the  aircraft  reaches  the  prescribed  altitude,  zero  dihedral 
and  maximum  pitch-up  elevator  angles  are  commanded.  These  signals  are  held  until  touch-down. 
Figure  20  shows  the  perching  signal  sent  at  the  0.6  s  mark.  The  angle  of  attack  builds  up  to 
30  deg,  causing  the  speed  to  reduce,  and  the  aircraft  climbs  momentarily.  Flight  speed  is  halved 
within  1  s  to  3  m/s.  After  a  brief  ascent,  the  MAV  lands  at  a  low  angle  of  attack.  It  is  interesting 
to  note  that  the  final  speed  has  reduced  substantially  even  without  using  wing  twist.  Addition  of 
wing  twist  would  not  only  enable  a  further  reduction  in  the  final  speed,  but  also  provide  for  better 
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roll  and  yaw  control  during  the  approach. 


(a)  a,  6  and  7  (b)  Flight  Speed 
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(c)  Trajectory 

Figure  20:  Flight  parameters  during  a  perching  attempt  that  was  triggered  at  1.5m  above  the 
ground 


9.8  Use  of  Trailing  Edge  Flaps  for  Mitigating  Control  Effectiveness  Problems 

Recall  that  the  sign  of  the  effectiveness  depends  on  the  sign  (xo,Cl/c  +  Cm^ac),  where  xa/c  is  the 
non-dimensional  distance  between  the  center  of  gravity  and  the  quarter-chord  line.  Furthermore, 
Cm,ac  <  0,  and  therefore,  at  small  angles  of  attack,  the  control  effectiveness  is  negative  and  it  is 
positive  at  higher  angles  of  attack.  For  an  intermediate  range  of  angles  of  attack,  the  sign  depends 
on  the  angular  rates  as  well. 

In  [26],  it  was  shown  that  one  way  to  get  around  this  problem  is  to  use  trailing  edge  (TE)  flaps. 
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TE  flap  deflection  leads  to  a  greater  increase  in  Cl  as  compared  to  the  reduction  in  CmAC.  It  is 
of  interest  to  find  the  flap  deflection,  as  a  function  of  a,  which  will  guarantee  a  certain  positive 
control  effectiveness.  Suppose  the  desired  value  is  0.025  (see  [26]  for  details).  Then,  for  the  aircraft 
considered  in  this  paper,  it  follows  that 

(~1 

+  CmtaC  +  0.72 5f  =  0.025,  (where  <5/  is  the  flap  deflection) 

0.07  + 0.5a -  0.1311  +  0.72,5/ =  0.025  =>  5f  =  0.12  -  0.69a  (40) 

The  benefit  of  a  uniformly  positive  control  effectiveness,  however,  comes  at  a  price:  the  aircraft 
is  forced  to  fly  in  a  high-lift,  high  drag  configuration  across  the  flight  envelope.  This  necessarily 
means  that  the  aircraft  will  fly  slower  than  normal.  However,  note  that  the  flight  path  angle  can  still 
be  controlled  effectively  using  symmetric  dihedral  deflection.  Experimental  results  which  validate 
this  concept,  and  demonstrate  its  effectiveness  in  a  successful  perching  maneuver,  are  reported  in 
[26].  Figure  21  shows  the  flight  parameters  measured  during  a  recently  conducted  set  of  perching 
experiments,  while  Fig.  22  shows  a  collection  of  snapshots  depicting  the  perching  maneuver. 


(a)  Trajectory  in  the  x  —  y  plane 


(b)  Trajectory  in  the  x  —  z  plane 
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Figure  21:  Flight  parameters  measured  during  a  recent  perching  experiment. 


37 


Figure  22:  Snapshots  of  a  perching  maneuver. 


10  PDE  Boundary  Control  of  Flexible  Wings 

The  motivation  for  considering  the  problem  of  boundary  control  of  beams  comes  from  the  problem 
of  controlling  flexible  wings  for  agile  aircraft  maneuvers.  Let  9(t,  y )  denote  the  twist  profile  of 
the  wing,  and  let  £(t,  y )  denote  the  bending  displacement  as  a  function  of  time  and  the  spanwise 
coordinate  y.  The  control  objective  is  to  ensure  that 


lim 

t—>  OO 


lim 

t—¥  OO 


lim 

t — yoo 


9(t,y)dy  =  0 


yO(t,y)dy  -  = 


(net  lift,  or) 

0  (net  rolling  moment,  and) 


£y(t,y)dy  -  R{t ) 


lim  (£(t,  L )  —  R(t))  =  0  (bending  displacement  of  the  tip), 

t~>  OO 


(41) 

(42) 

(43) 


where  H(t),  Hi(t)  and  R(t)  are  sufficiently  smooth  time- varying  signals.  Although  we  state  asymp¬ 
totic  convergence  as  the  objective,  we  will  prove  exponential  convergence  to  a  uniform  ultimate 
bound.  The  term  £y(t,y)dy  =  £(L)  is  a  measure  of  the  effective  wing  dihedral  which  is  a  key 
yaw  control  parameter  [62],  and  it  measures  the  amount  of  side  (y-)  force  produced  by  the  wing 
which,  in  turn,  produces  a  yawing  moment  on  the  aircraft. 

Remark:  An  important  question  that  concerns  systems  described  by  PDEs  is  that  of  well- 
posedness.  Well-posedness  of  the  closed  loop  systems  considered  here  can  be  shown  by  proving 
that  the  input-output  map  of  the  system  is  bounded  [82,  11],  For  the  twisting  dynamics  actuated 
by  root  control,  this  is  achieved  by  designing  the  control  to  map  the  system  onto  well-posed  and 
exponentially  stable  dynamics.  For  twisting  dynamics  actuated  at  the  wing  tip,  the  input-output 
map  is  actually  a  finite  order  ODE,  and  its  well-posedness  follows  from  the  standard  existence  and 
uniqueness  theorems  for  ODEs.  Finally,  the  well-posedness  of  the  closed  loop  bending  dynamics 
can  be  shown  using  a  transfer  function  approach  [11], 


38 


10.1  Boundary  Control  of  Twisting  Motion:  Root  Control 

PDE  backstepping  is  used  for  designing  a  boundary  control  law  for  wing  twisting  dynamics. 

Backstepping  is  carried  out  in  two  steps:  (a)  the  target  dynamics  are  identified,  and  (b)  a  back- 
stepping  (Volterra)  transformation  converts  the  system  dynamics  (in  this  case,  the  error  dynamics) 
into  the  target  dynamics  and  the  control  signal  u(t )  is  obtained  in  the  process.  The  method  de¬ 
scribed  by  Krstic  and  Smyshlyaev  [42]  is  used  here.  Let  9  denote  the  error  between  the  system 
state  and  the  desired  steady  state  value,  i.e.,  9  =  9  —  6^. 

Next,  consider  the  target  dynamics  described  by  the  PDE 

vtt  ~  bvtyy  -  avyy  =  (M  -  ap)v  -  bpvt, 

v(t,  0)  =  vy(t,  L)  =  0  (44) 

Using  the  method  of  separation  of  variables,  it  is  easy  to  check  that  the  eigenvalues  of  this  system 
are  the  solutions  of 


A2  +  b(v2  +  p)  A  +  (a(^2  +  p)  —  M)  =  0 
2n  +  1  7T 


where  v  = 


n  =  0,  ,1,  2,... 


2  L 

The  target  dynamics  are  stable  if  and  only  if  the  control  design  parameter,  p,  satisfies 

f  M  TV2  IT2  \ 

p  >  max  (  —  4^2’  -U2) 


(45) 


(46) 


A  dummy  spatial  variable  x  is  introduced  and  the  Volterra  transformation  between  9  (the  error 
between  the  actual  dynamics  and  the  desired  steady  state  in  equation  (44))  and  v  (the  target 
dynamics  for  the  error  state)  is  given  by 


v(t,  y )  =  9(t ,  V)  ~  J  x)9(t,  x)dx  (47) 

It  is  helpful  to  recall  that  the  9  dynamics  are 

9tt  -  b9tyy  -  a9yy  =  MO,  0 (t ,  0)  =  u(t),  6y(t,  L)  =  0  (48) 

In  order  to  solve  for  k(x,y),  substitute  equations  (47)  and  (48)  into  equation  (44).  Next,  isolating 
the  coefficients  of  v  and  vt,  the  following  Klein  -  Gordon  PDE  for  k(y,x)  is  derived  [42] 


kxx(y,x)  -  kyy(y,x)  =  —pk(y,  x) 
k{y,y)  =  |(L  -  y),  kx(y,L)  =  0 

The  control  input  is  found  from  equation  (47) 


u(t)  =  0(t,  0)  =  —  f  k(0,x)9(t,x)dx 

Jo 


(49) 


(50) 
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The  solution  is  given  in  terms  of  the  modified  Bessel  function  I\  *: 


k(y,x)  =  p(L  -  y ) 


HVm  -  y )2  -{L-  a)2)) 
Vp((L~y)2  -  (L~x)2)) 


(51) 


A  few  observations  are  worth  noting  here. 


1.  If  the  wing  is  reasonably  stiff 


Mir, 


M 


7T' 


A 


GJ  a  ML2)' 
be  stabilized  with  p  =  0,  i.e. ,  with  no  additional  control  input. 


=  —  <  — g  ;  i.e.,  GJ  >  4:L2MIp/ir  ,  the  system  can 


2.  For  stability,  it  is  essential  that  b  >  0,  i.e.,  the  Kelvin- Voigt  damping  coefficient  is  always 
positive.  A  negative  damping  could  be  introduced  due  to  aerodynamics,  but  it  can  be  com¬ 
pensated  by  the  term  bpvt  and  wing  flutter  can  be  prevented.  The  compensation  in  damping 
imposes  an  additional  constraint  on  p. 


3.  The  controller  in  equation  (50)  requires  that  the  twist  angle  at  all  points  on  the  wing  be  known. 
This  difficulty  can  be  circumvented  by  designing  a  PDE-based  observer  [42]  or,  practically, 
using  a  series  of  distributed  sensors  and  fitting  their  output  with  an  a  priori  designed  spline. 


4.  Damping  and  stiffness  cannot  be  added  independently.  They  are  added  in  the  ratio  b/a. 

5.  Finally,  the  control  law  does  not  require  that  a,  b  or  M  be  known  for  the  purpose  of  regulation. 
We  only  need  to  know  bounds  on  a  and  M  to  choose  an  appropriate  value  of  the  gain  p. 


The  backstepping  controller  shown  above  is,  strictly  speaking,  a  stabilizing  controller:  it  regu¬ 
lates  6(t,y )  —>  0.  The  problem  of  designing  a  tracking  controller  which  achieves  the  objective  in 
(41)  needs  further  work.  Consider  the  system  realized  via  the  coordinate  transformation 

/y 

6(t,x)dx  (52) 

Physically,  w(t,  y )  measures  the  lift  generated  by  the  outboard  section  of  the  wing  starting  at  y  and 
terminating  at  the  wing  tip.  The  lift  produced  at  the  wing  tip  is  zero,  which  is  physically  correct. 
Note  that  0y(t,L )  =  0  at  the  free  end  y  =  L.  Hence,  it  follows  that 

ry  rv 

wtt(t,  y)  =  /  0t.t(t,x)dx=  /  (bOtxx{t,  x)  +  aOxx(t,x)  +  M9(t,x))dx 
Jl  Jl 

=  bOty(t,  y)  +  a0y{t ,  y)  +  Mw(t ,  y)  =  bwtyy(t,  y)  +  awyy(t ,  y)  +  Mw(t,  y)  (53) 


Thus,  the  dynamics  of  w  are  described  by  the  PDE 


wu  —  bwtyy  —  awyy  =  Mw,  w(t,  L)  =  0,  wy(t,  0)  =  u(t),  (54) 

where  u(t)  is  the  control  signal  (Or),  and  recall  that  a  =  GJ/Ip  and  b  =  pa.  Furthermore,  we  defined 
M  so  that  M0  =  —xacFb/ Ip,  where  is  a  linear  function  of  6.  Note  that  wy(t,y )  =  6(t,y).  Recall 
that  the  control  objective  is  to  ensure  that  limt_>.oo(J0  Ody  —  H(t))  =  0  (see  equation  (41)).  The 

*A  modified  Bessel  function,  In(y),  satisfies  y2I'n(y)  +  yl'n{y)  —  ( y 2  +  n2)In(y)  =  0. 
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Figure  23:  Block  diagram  showing  the  tracking  controller  for  twist  (8)  dynamics.  This  structure  is 
identical  to  the  classic  strict  feedback  structure  for  systems  described  by  ODEs  [41]. 


control  objective  is  now  recast  to  ensuring  that  limi^.00(t(;(t,  0)  +  H(t ))  =  0,  i.e.,  lim^oo  w(t,  0)  = 

-m- 

The  method  for  designing  a  tracking  controller  involves  the  following  three  steps: 

1.  Obtain  a  backstepping  transformation  w  >— >  v,  although  the  boundary  conditions  are  changed 
to  match  the  tracking  requirement, 

2.  Identify  the  boundary  conditions  for  tracking  and 

3.  Derive  a  motion  planning-based  design  for  the  boundary  conditions  of  the  v  dynamics. 

The  details  of  this  method  can  be  found  in  [63].  The  control  design  procedure  has  been  illustrated 
via  a  block  diagram  in  Fig.  23. 

10.2  Boundary  Control  of  Twisting  Motion:  Wing  Tip  Control 

As  in  the  previous  section,  one  can  design  a  backstepping  controller  for  the  case  where  a  control 
moment  is  applied  to  the  free  end  (y  =  L )  of  the  wing  while  the  other  end  (y  =  0)  is  clamped.  In 
fact,  the  procedure  in  both  cases  is  identical,  although  the  final  expressions  for  the  control  law  differ 
slightly.  Alternately,  in  case  of  MAVs,  one  may  do  without  a  stabilizing  controller.  The  “tracking 
half”  of  the  controller  (without  the  stabilizing  backstepping  controller)  may  be  designed  using  the 
output  measurements.  This  method  is  useful  for  adaptive  designs  as  well.  We  consider  the  wing 
model 

Ott  ~  b8tyy  -  aOyy  =  M0 ,  0 (t ,  0)  =  0,  0y(t,  L)  =  u(t )  (55) 

where  the  control  input  is  a  moment  applied  at  the  wing  tip.  The  control  objective  is  to  ensure 
that 

hm  (  f 

t-xx>  yj q 


9(t,y)dy  -  H(t))  =  0 
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Let  e(t)  =  /0L  0(t,y)dy  —  H(t )  denote  the  error  which  needs  to  be  regulated.  Then, 

e=  f  Oa(t,y)dy  -  H(t) 

Jo 

L 

(aOyy  +  bdtyy  +  MO)  dy  -  H{t ) 

=  aOy{L)  -  aOy{ 0)  +  bOty(L)  -  b0ty( 0)  +  Me(t)  -  H(t )  +  MH(t) 

=  bu(t)  +  au(t)  —  a0y(  0)  —  b0ty(  0)  +  Me(t)  —  H(t)  +  MH(t )  (57) 

A  dynamic  controller  of  the  form 

6ti(f)  +  art(t)  =  H(t)  —  MH(t )  —  (M  +  k)e(t)  —  kce(t)  +  a0y(O)  +  60jy(O)  (58) 

renders  the  system  into  the  spring-mass  form 

e(i)  +  kce(t )  +  fce(t)  =  0.  (59) 

The  control  law  in  equation  (58)  suggests  that  0  need  not  be  monitored  or  measured  at  all  locations 
on  the  wing.  Instead,  only  0y( 0)  needs  to  be  measured  or  estimated.  The  reference  signal  H(t )  is 
known.  It  may  be  difficult  to  inject  damping  because  e(t)  is  the  rate  of  change  of  the  lift  and  in 
practice,  would  require  differentiating  noisy  acceleration  signals. 

Another  interesting  observation  is  that  although  the  PDE  system  had  an  infinite  relative  degree 
when  the  root  twist  was  chosen  as  the  control  input,  the  relative  degree  is  2  when  twisting  moment 
at  the  wing  tip  is  considered  as  the  input.  This  facilitates  the  control  law  design  in  this  section 
considerably.  The  control  law  design  described  is  this  section  lends  itself  readily  to  adaptation 
should  a  and/or  M  be  unknown. 

Consider  the  dynamics  in  equation  (55)  with  the  objective  in  equation  (56)  and  suppose  that 
M  is  unknown.  The  control  law  in  (58)  is  modified  so  that 

bu(t)  +  au(t )  =  H(t)  —  M(t)(H(t)  +  e(t))  —  ke(t)  —  kce(t )  +  a0y(  0)  +  b0ty(  0)  (60) 

where  M(t )  is  the  estimated  value  of  M.  The  error  dynamics  are  described  by  the  ODE 

e(t)  +  kce[t )  +  ke{t)  =  +  e(t)),  (61) 

where  M{t )  =  M(t)  —  M(t).  The  adaptive  law  for  M(t) 

M(t )  =  7  Proj  ^M(f),  xTP  J  (e(f)  +  H)^j  (62) 

ensures  that  the  error  e(t)  remains  bounded  on  the  order  O(lf'y),  so  that  the  adaptive  gain  can  be 
chosen  to  match  the  tracking  accuracy  requirements.  Although  a  and  b  were  assumed  to  be  known, 
the  aforementioned  analysis  can  be  repeated  to  accommodate  an  unknown  a  and  b  as  well. 
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10.2.1  Tracking  and  Stability 


The  problems  of  tracking  and  stabilization  are  distinct  because  the  PDE  system  is  infinite  dimen¬ 
sional.  Nevertheless,  a  tracking  controller  improves  stability  as  described  presently.  The  best  way 
to  understand  the  influence  of  a  tracking  controller  is  to  set  H  =  0.  Make  a  coordinate  transforma¬ 
tion  w(t,  y )  =  [/  6(t,  y)dy  so  that  achieving  H  =  0  is  equivalent  to  achieving  w(0)  =  w(t,  0)  =  0.  It 
follows  that  the  transformed  dynamics,  with  H  =  0  ensured  by  the  tracking  algorithm,  are  given 
by 

UHt(t,y)-bwtyy(t,y)-aWyy(t,y)  =  Mw(t,y)-bwtyy(t,0)-awtyy(t,0),  w(L)  =  w( 0)  =  0,  u/(0)  =  0 

(63) 

The  third  boundary  condition  is  not  entirely  independent.  Let  w(t,y)  =  r](t)(j)(y).  Then,  we  get 

ij(t)  -  My(t)  =  4>"(y)  -  r(0)  = _ \  2 

br,(t)  +  ar,(t)  <t>{y)  ’  1  j 

where  A  is  a  constant.  It  follows  that  the  condition  for  stability  is  M/a  <  A2.  The  differential 
equation  for  </>{y)  can  be  solved  to  get 

(j){y)  =  A  sin(A y)  +  B  cos(A y)  +  </>"(0)/A2  (65) 

The  boundary  conditions  </>(0)  =  4>'(0)  =  0  lead  to  4>(y)  =  B(cos(Xy)  —  1).  Finally, 

<j>(L)  =  0  =>  A  =  2nn/L,  (j>(L)  =  0,  ra  =  l,2, ...  (66) 

Had  we  not  imposed  the  condition  H( 0)  =  0,  we  would  have  obtained  M/a  <  7t2/(4L2)  as  the 
condition  for  stability.  Since  the  condition  for  open  loop  stability  is  given  by  M/a  <  A2  =  4tt2/L2, 
it  follows  that  the  stability  margin  improves  by  a  factor  of  sixteen  using  only  the  tracking  controller, 
although,  it  does  not  stabilize  the  wing  for  all  values  of  M  and  a  as  backstepping  does.  In  principle, 
the  tracking  controller  converts  the  wing  from  a  cantilever  beam  to  a  clamped-clamped  beam.  In 
practice,  this  translates  to  the  ability  to  increase  the  wing  flexibility  by  an  order  of  magnitude,  or 
increase  the  wing  divergence  speed  four-fold.  This  is  shown  in  Fig.  24. 


10.3  Control  of  Rolling  Moment 

An  abstract  measure  of  the  rolling  moment  is  f  yO(t,  y)dy,  defined  in  equation  (42).  Let  e/(t,  y)  = 
/0L  y9[t,  y)dy  —  H[(t ),  where  Hi(t)  denotes  the  reference  value  of  the  rolling  moment  to  be  tracked, 
and  ei(t)  is  the  tracking  error.  Differentiate  e/(t)  twice  with  respect  to  time: 


ei(t)=  ydu(t,y)dy- Hi(t) 


nL  rL 

=  /  y{bdtyy(t,y)  +a0yy(t,y))dy  +  M  \  y9(t,y)dy  -  Ht(t) 


(67) 


/  o  Jo 

=  L(bOty(t ,  L)  +  aOy(t,  L))  -  b(et(t ,  L)  -  9t(t,  0))  -  a(0(t,  L)  -  6{t ,  0))  +  Mej(t)  +  MHt(t)  -  Hi(t ) 


An  interesting  observation  is  that  ei(t )  has  a  relative  degree  of  2  with  respect  to  0y(L)  (tip  control) 
as  well  as  0(0)  (root  control).  In  particular,  it  means  that  a  considerably  simpler  controller  than 
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Figure  24:  Torsional  divergence  speed  as  a  function  of  the  Young’s  modulus  of  the  wing.  The  blue 
curve  is  the  open  loop  divergence  speed,  while  the  red  curve  is  the  divergence  speed  after  adding 
the  tip-based  twist  controller. 


the  backstepping  controller,  on  the  lines  of  the  adaptive  controller  in  the  previous  section,  can  be 
implemented  for  a  root-based  actuator.  Indeed,  in  aircraft,  the  lift  is  controlled  using  the  horizontal 
tail,  while  wing-based  flaps  (ailerons  and  spoilers)  are  used  primarily  for  roll  control. 

10.4  Root  Control  for  Wing  Bending 

The  bending  PDE  can  be  written  as 

Fl^tt  T  VEIb^tyyyy  T  F Ib^yyyy  =  En(t,  @tt)i 

£(*,0)  =  Zyy{t,L)  =  €yyy(t,L)  =  0,  £y (t,  0)  =  u(t) ,  (68) 

where  the  control  input  u(t)  =  6n(t)  is  designed  to  ensure  that  £(i,  L)  =  R(t)  in  (43).  The 
acceleration  term  corresponding  to  mxec6  (due  to  twist)  has  been  moved  to  the  right  hand  side 
and  merged  into  F^,  so  that  Fn  =  Ff,  +  mxec6.  The  interesting  point  about  equation  (68)  is  that 
the  right  hand  side  is  independent  of  £,  and  6  is  an  average  value  obtained  from  the  faster  twist 
dynamics.  Therefore,  unlike  the  twisting  dynamics,  the  onset  of  instability  in  the  bending  dynamics 
will  correspond  to  the  damping  becoming  negative. 

Bending  dynamics  are  usually  stable,  in  that  there  is  no  direct  aerodynamic  driver  of  instability 
as  there  is  in  the  case  of  torsion.  If  the  beam  dynamics  are  unstable  PDE  backstepping  can  be  used 
to  inject  damping  into  the  system.  A  procedure  for  backstepping,  based  on  a  transformation  of  the 
bending  dynamics  into  the  Schrodinger  equation,  has  been  detailed  in  Krstic  and  Smyshlyaev  (See 
Ch.  8  of  Ref.  [42]).  A  tracking  controller  can  be  added  on  top  of  the  backstepping  controller  for 
tracking.  We  focus  on  a  motion  planning-based  design  for  a  tracking  controller.  However,  note  that 
F(t,  y )  is  usually  not  measurable  in  practice,  although  one  may  estimate  its  spatial  profile  from  the 
wing  geometry  [49]. 

In  this  section,  we  design  an  observer-based  controller  to  facilitate  a  motion-planning-based 
tracking  controller  for  bending.  The  perturbation  observer  does  not  predict  the  system  states.  It 
uses  projection-based  adaptation  to  estimate  F(t,  y )  which  would  be  unknown  in  practical  situa¬ 
tions.  The  perturbation  observer  is  split  into  a  “particular”  and  a  “homogeneous”  component  (the 
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notions  will  be  made  more  precise  in  this  section).  In  particular,  the  homogeneous  component  is 
stable  and  not  driven  directly  by  external  feedback.  Thus,  it  is  simpler  to  design  a  control  law  for 
it.  The  same  control  signal  is  sent  to  the  actual  system,  whose  states  then  converge  exponentially 
to  a  bounded  envelope  around  the  observer  states. 

Let  F(t,  y)  =  W (t)T (j>b(y)  +  &(t)  where  W(t)  and  a(t)  are  unknown  and  bounded  with  known 
bounds.  The  set  of  functions,  4>b(y),  can  be  chosen  to  get  a  satisfactory  bound  on  a,  and  using  a 
knowledge  of  the  wing  geometry  [49]. 

The  perturbation  observer  can  be  designed  as  a  sum  of  two  states,  £  =  £/,  +  where  the 
dynamics  of  the  two  states  £/,_  and  £p  are  described  by  the  following  PDEs: 

U )  T  bbCpttyyyy(t,  y)  +  Clbip,yyyy(tj  y)  =  y)  dbP£,p{t j  y)  +  W ( t )  +  O'(f), 

£ p,yy( L)  =  ^p,yyy{F)  =  £p(0)  =  £p,?;(0)  =  0  (69) 

for  the  particular  component,  where  £p  =  Cp  ~  £>  and 

£ h,tt  T  bb^h,tyyyy  T  (!h^h:yyyy  —  bbP^h,t  f,hP/h  i 

ih,yy(L )  =  £h,yyy(L)  =  £h(0)  =  0,  £h,y{ 0)  =  u(t)  (70) 

for  the  homogeneous  component.  Here  ^  ^  — ^  and  p  >  0  is  chosen  to  ensure  desirable  convergence 

properties.  Projection-based  adaptive  law  for  W(t)  and  a (t), 

W(t)  =  7aProj  (w(t),  -J  (it  +  Si)4>b(y)dy 

°(t)  =  7«Proj  (it  +  Si)dy\  (71) 

lead  to  error  dynamics  between  the  perturbation  observer  and  the  actual  system  that  are  stable, 
and  in  fact,  globally  uniformly  bounded. 

This  observer  can  be  used  in  the  motion  planning  algorithm  to  ensure  that  the  wing  tip  tracks 
the  desired  displacement  profile.  The  motion  planning  is  performed  in  practice  on  the  homogeneous 
half  of  the  observer.  Note  that  the  homogeneous  half  is  a  stable  system  and  its  dynamics  are  not 
influenced  in  any  way  by  £  or  £p.  The  details  of  motion  planning  are  given  in  [63].  Figure  25  is  a 
block  diagram  of  the  perturbation  observer.  Notice  that  the  reference  input,  R(t)  from  (43),  enters 
only  the  homogeneous  component  of  the  perturbation  observer,  while  the  feedback  from  the  actual 
system  only  enters  the  particular  component. 

10.5  Simulations 

Simulations  are  carried  out  in  Matlab  using  a  Galerkin-based  approach  to  convert  the  PDE  system 
into  ODEs.  The  Galerkin  truncation  is  not  used  as  a  basis  for  control  law  design,  so  no  danger 
of  a  “spillover  instability”  arises.  Figure  26  demonstrates  the  regulation  of  twist  dynamics  using 
the  backstepping  controller  derived  in  equation  (50),  with  the  transformation  in  equations  (47)  and 
(49).  The  value  of  M/a  was  set  to  8,  where  a  =  GJ/Ip.  A  value  of  p  =  4  yielded  an  unstable 
response,  while  the  response  was  stable  for  p  =  8.  Recall  the  following  condition  for  stability  with 
L  =  1:  p  >  M/a  —  7r2/ 4  ~  2.53.  The  backstepping  controller  works  even  when  M (y)  =  M ( 1  —  y2) 
is  used  (to  mimic  an  elliptical  lift  distribution  over  the  wing)  instead  of  a  constant  M(y)  =  M. 
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Figure  25:  Block  diagram  for  the  perturbation  observer  coupled  to  the  system  dynamics.  The 
control  signal,  u(t),  is  generated  from  motion  planning,  while  R(t)  is  the  desired  reference  signal 
from  (43). 


(a)  Unstable  response  with  p  =  4  (b)  Stable  response  with  p  —  8  for  spatially-varying 

M(y) 


Figure  26:  Regulation  of  the  twist  dynamics  using  the  backstepping  controller  in  equation  (50), 
with  the  transformation  in  equations  (47)  and  (49).  The  plots  were  obtained  for  M/a  =  8,  while  p 
was  increased  to  ensure  stability.  Each  plot  is  a  collection  of  snapshots,  where  the  lines  get  darker 
with  time. 
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The  backstepping  controller  can  be  added  on  top  of  a  tracking  controller. 

Figure  27  shows  the  simulation  of  a  wing  actuated  by  tip  control.  The  value  of  M/a  was  chosen 
so  that  stability  is  assured  without  the  need  for  a  dedicated  stabilizing  controller.  The  first  plot  was 
obtained  for  a  for  a  system  where  the  aerodynamics  were  assumed  to  be  linear  but  unknown.  The 
second  plot  was  obtained  for  the  case  where  the  aerodynamics  were  additionally  spatially  varying 
and  rendered  the  open  loop  dynamics  unstable.  The  time  histories  of  the  tracking  error  e(t),  and 
the  control  signal  u(t),  for  the  second  case  are  plotted  in  Fig.  28.  In  both  cases,  the  twist  amplitude 
converges  to  the  steady  state  value  with  satifactory  transients.  The  error  metric  equation  (41)  is 
also  seen  to  be  very  small.  Note  that  the  control  signal,  obtained  from  Eq.  (60),  is  almost  noise- 
free.  This  is  because  the  dynamic  controller  acts  like  a  low  pass  filter  and  ensures  that  noisy  terms 
arising  from  the  high  gain  terms  on  the  right  hand  side  of  Eq.  (60)  are  filtered  out. 

Finally,  figure  fig:BendSims  shows  the  time  histories  tip  displacement  for  a  wing  whose  bending 
motion  is  actuated  by  a  root-based  actuator.  The  dashed  lines  show  the  reference  signal  for  wing 
tip  displacement,  while  blue  lines  show  the  actual  displacement  of  the  wing  tip.  The  first  simulation 
was  obtained  for  the  case  where  F(t,  y )  was  set  to  zero,  while  the  second  plot  was  obtained  for 
a  time  varying  F(t,y).  The  motion  planning  algorithm  used  a  seventh  order  polynomial  in  the 
spatial  variable  y.  In  both  cases,  the  tracking  is  seen  to  be  satisfactory  and  noise  free. 


(a)  M  unknown,  constant,  open  loop  stable  (b)  Spatially  varying,  unknown  M(y),  unstable  dy¬ 

namics 

Figure  27:  Twist  profile  of  the  wing  as  a  function  of  time  when  the  the  adaptive  controller  in 
equation  (60)  is  applied  at  the  wing  tip.  Three  cases  have  been  examined  here,  with  9(t ,  y)dy  = 
0.05  as  the  desired  output.  Each  plot  shows  appropriately  chosen  snapshots,  with  lines  getting 
darker  with  time. 


10.6  Section  Summary 

Tools  from  PDE  boundary  control  show  considerable  promise  for  the  efficient  control  of  flexible 
structures  in  general,  and  of  flexible  wings  in  particular.  The  key  contributions  of  this  work 
were:  (a)  the  design  of  a  PDE  backstepping  based  tracking  controller  for  wing  twist,  (b)  the 
finite  degree  problem  formulation  for  wing  twist  control  using  a  tip-based  actuator,  (c)  the  study 
of  the  enhancement  in  stability  due  to  a  tracking  controller,  and  (d)  the  design  of  a  perturbation 
observer  based  controller  for  wing  bending.  It  is  safe  to  claim  that  every  contribution  listed  above 
leads  to  a  related  open  problem  in  PDE  control.  The  most  important  open  problem  is  the  design 
of  a  coupled  twist-bending  controller. 
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(a)  Time  history  of  e(t)  (b)  Time  history  of  u(t) 

Figure  28:  Time  histories  of  the  error  e(t)  and  the  control  signal  u(t )  for  the  case  of  spatially 
varying,  unknown  M(y )  with  unstable  open  loop  dynamics. 


Figure  29:  Time  history  of  the  wing  tip  displacement  for  two  classes  of  reference  signals:  a  pulse 
and  a  sine  wave.  The  first  two  plots  were  obtained  with  the  right  hand  side  set  to  zero,  while  the 
third  plot  was  obtained  for  a  time  varying  F(t,y)  in  equation  (68). 

11  Conclusions 

The  objective  of  this  work  was  to  study  the  unique  dynamics  and  control  characteristics  of  tailless 
MAV  equipped  with  articulated  wings.  We  studied  the  performance  and  stability  using  a  combi¬ 
nation  of  literal  approximations,  numerical  trim  and  stability  computations,  and  experiments.  We 
showed  that  the  dihedral  angle  of  the  wing  can  be  varied  symmetrically  to  obtain  an  additional 
degree  of  freedom,  namely  the  ability  to  change  flight  path  angle  independently  of  the  flight  speed. 
We  demonstrated  that  asymmetric  dihedral  settings  can  be  used  to  perform  rapid  turns  and  control 
the  sideslip.  From  the  standpoint  of  control,  the  most  important  observation  was  the  discovery  of 
maneuver-dependent  control  effectiveness  reversal. 

Bat-like  flight  is  a  challenging  problem  that  cannot  be  solved  via  averaging  or  with  traditional 
tail-derived  stability.  We  have  demonstrated  the  ability  to  stabilize  and  control  longitudinal  motions 
via  CPGs  with  the  RoboBat.  As  expected,  the  top-level  controllers  are  of  low  dimension  and  can 
be  made  very  simple,  because  most  of  the  hard  work  is  done  by  the  CPGs.  Given  the  extent  of 
mechanical  coupling  in  the  design,  it  is  remarkable  that  such  control  was  immediately  as  effective 
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as  it  was.  Further  work  can  still  be  done  to  create  a  pattern  generator  layer  so  to  optimize  the 
output  waveforms.  Additionally,  we  expect  to  quantify  the  forces  and  moments  actually  produced 
via  the  aerodynamic  model,  so  that  we  can  make  better  predictions  for  a  free-flying  robotic  bat. 

As  mechanical  design  of  actuators  develops,  we  expect  robotic  fliers  in  free  flight  to  be  able  to 
utilize  the  key  feature  of  phase  synchronization  and  control  of  phase  differences  in  stability  and 
control  of  body  motions.  The  major  problem  of  identifying  a  method  of  proving  such  stability 
analytically  is  still  open.  However,  this  work  has  demonstrated  the  successful  CPG-based  flight 
control  result  experimentally  by  using  the  RoboBat.  Since  this  CPG  controller  design  also  features 
rapid  inhibition  of  oscillation,  it  leads  to  the  important  problem  of  gliding  flight  and  maneuvers 
while  gliding. 

The  performance  and  stability  characteristics  of  a  flexible  aircraft  were  compared  with  those 
of  a  rigid  aircraft,  assuming  that  the  wings  of  the  flexible  aircraft  were  statically  deformed.  We 
presented  a  metric  called  the  effective  dihedral  for  flexible  wings  which  allows  the  results  from  the 
analysis  of  a  rigid  aircraft  to  be  extended  to  flexible  aircraft.  Moreover,  it  allows  us  to  identify 
the  extent  to  which  conclusions  regarding  the  performance  of  rigid  aircraft  apply  to  a  flexible¬ 
winged  aircraft.  Although  we  did  not  observe  any  difference  in  stability  characteristics,  there  were 
interesting  differences  in  the  turn  performance.  Although  wing  flexibility  was  shown  to  help  reduce 
the  sideslip  significantly,  it  was  also  shown  to  reduce  the  maximum  attainable  turn  rate  under 
a  zero-sideslip  constraint  compared  to  a  rigid  aircraft  with  identical  elevator  settings.  This  was 
attributed  to  the  reduction  in  flight  speed  due  to  wing  twist. 

The  open  loop  characteristics  of  a  tailless  MAV  were  tested  experimentally  to  verify  the  ana¬ 
lytical  predictions.  Control  laws  motivated  by  dynamic  inversion  were  designed  and  tested  on  an 
experimental  aircraft.  The  experiments  exposed  the  impediments  that  arise  due  to  the  reversal  of 
control  effectiveness.  At  the  same  time,  they  demonstrated  the  ability  of  articulated  wings  to  aid 
steep  descents  and  perching  maneuvers,  as  well  as  their  capability  for  yaw  control. 

Finally,  we  designed  PDE-based  control  laws  for  controlling  the  deformation  of  a  flexible  wing 
to  achieve  a  net  aerodynamic  force  or  moment.  We  considered  cases  where  the  actuators  were 
based  at  the  wing  root  as  well  as  the  tip.  PDE  backstepping-based  control  laws  were  developed 
for  controlling  wing  twist  using  root-based  actuators.  We  showed  that  a  tracking  controller  could 
bring  about  a  significant  improvement  in  the  stability  margin  of  the  wing  dynamics,  measured 
by  the  critical  flight  speed  or  elastic  modulus  for  the  onset  of  instability.  We  showed  that  the 
time  scale  separation  between  the  bending  and  twisting  dynamics  depends  primarily  on  the  aspect 
ratio  of  the  wing  and  the  flight  speed.  We  designed  a  controller  for  bending  independently  of 
the  one  for  wing  twist.  The  controller  designed  for  bending  used  a  novel  idea  based  on  splitting  a 
perturbation-observer  into  two  parts,  one  of  which  accommodated  the  external  forces  and  the  other 
which  accommodated  the  boundary  control.  Thereafter,  motion  planning  was  used  to  design  the 
boundary  controller,  with  the  understanding  that  it  could  be  added  readily  on  top  of  a  stabilizing 
controller  if  required. 
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